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We have developed a semiclassical approach to solving the Bogoliubov - de 
Gennes equations for superconductors. It is based on the study of classical orbits 
governed by an effective Hamiltonian corresponding to the quasiparticles in the 
superconducting state and includes an account of the Bohr-Sommerfeld quanti- 
sation rule, the Maslov index, torus quantisation, topological phases arising from 
lines of phase singularities (vortices), and semiclassical wave functions for multi- 
dimensional systems. The method is illustrated by studying the problem of an 
SNS junction and a single vortex. 



1. INTRODUCTION 

A most convenient microscopic theory of superconductivity is provided 
by the Bogoliubov-de Gennes (BdG) equations H. On the one hand it 
encapsulates all the basic ideas of Bardeen, Cooper and Schrieffer and on 
the other hand it is the general form of the Kohn-Sham Euler-Lagarange 
equations of the density functional theory for superconductors. Ever since 
its formulation semiclassical methods have been used, very successfully, for 
investigating its solutions under various circumstances pi pi 0]. In this 
paper we wish to contribute to making this method even more effective by 
developing further the particular approach, based upon effective classical 
orbits, pioneered by Azbel' j|. 

As in the case of the standard Schrodinger equation there are two semi- 
classical approaches to the problem of solving the BdG equations. The first, 
initiated by Andreev j^, |j, de Gennes Q, and Bardeen and co workers 
||, utilises the WKB form for the wave function and takes advantage of 
the slowly varying amplitudes to convert the second order equation into 
one that is first order in the derivatives and hence more readily soluble. 
However the usual machinations of Schodinger quantum mechanics e.g. of 
matching logarithmic derivatives etc. are still employed. 
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The second approach, initiated by Azbel' Q, uses effective classical 
Hamiltonians and orbits together with Bohr quantisation conditions to 
study the quasiparticle energy spectra and the corresponding wave func- 
tions. We extend Azbel's approach by including an account of complex 
order parameters, the origin of the Maslov indices, the construction of a 
3-dimensional wave function, and using torus quantisation, that is to say 
the full machinery of modern semiclassics ||, ]j| [| . 

The principle motivation behind such systematic development of semi- 
classical methods for superconductors is the need to adopt this powerful 
technique to the treatment of Type II superconductors and those with 
exotic, p- and d-wave, pairing. This need arises from the difficulty of solv- 
ing numerically differential equations, such as the BdG equations, which 
feature many wildly different length scales, such as lattice parameters, a, 
coherence length, £, penetration depth, radii of Landau orbits, and flux 
lattice unit cell sizes. We hope that the replacement of the numerical prob- 
lems of integrating differential equations by solving Hamiltons equations 
for classical orbits can alleviate some of these difficulties. 

The layout of this paper is then as follows. In section || we present a 
multicomponent semiclassical theory for the BdG equations. This starts 
with a discussion of the general form of the semiclassical ansatz before 
showing how the BdG equations can be reduced, at zeroth order, to the 
solution of a pair of Hamiltonian systems describing the dynamics of quasi- 
particles along orbits comprised of both particle-like and hole- like segments. 
We also include a discussion of topological phases arising from singularities 
in the phase of pairing potential. We then solve the first order equations 
and construct the general multicomponent wave function. We briefly dis- 
cuss torus quantisation and demonstrate that EBK quantisation conditions 
(a generalisation of Bohr-Sommerfeld quantisation rules) can only be con- 
structed in the absence of the above mentioned singularities. In section ^| 
we construct an effective semiclassical theory whose Hamiltonian systems 
are ^.-dependent. By extending the semiclassical theory in this way we 
remove the obsticles to constructing the EBK rules and consequently we 
present a generalised EBK rule which includes the contribution from both 
the Maslov index and topological phases. We conclude the section with a 
discussion of the interpretation of a theory based upon ^.-dependent Hamil- 
tonian systems. The remaining parts of the paper are devoted to two appli- 
cations. Thus in section ^ we apply the theory to a superconductor-normal 
metal-superconductor junction and include a discussion of Andreev retro- 
reflection and the semiclassical spectrum. In an accompanying appendix 
we calculate the Maslov index by analytically continuing the semiclassical 
wave functions taking into account of Stokes phenomenon. Then in section 
U we investigate a single vortex which is an ideal example for demonstrating 
the role of the ^-dependence in the Hamiltonians and the consequences of 
the phase singularities in the pairing potential. The paper concludes with 
a discussion and summary. 
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2. SEMICLASSICAL THEORY FOR SUPERCONDUCTORS 



2.1. General form for the multicomponent WKB ansatz 

To begin, let us recall the conventional semiclassical (2[ ||, [j], [j| BCS 
theory of the superconducting state. The context for the most general 
version of this is provided by the Bogoliubov-de Gennes equations for the 
two component wave function for a quasiparticle: 

H(p,r) |A(r)|e^ r ) \ / u A (r) \ _ / u x (r) \ 
|A(r)|e-^W -H*(p,r) J { u A (r) ) ~ * x I v x {v) )> (ZA) 



where 



H(p, r) = -L (p + eA(r)) 2 + U(r) - e F , (2.2) 
2m 



6f is the Fermi energy, U(r) is an external potential and A(r) = | A(r)|e l ^^ r -' 
is the complex pairing potential satisfying the self-consistency condition 

A(r) = gJ2 ux(r)v* x (r) (1 - 2f{E x )) . (2.3) 

A 

(Here g is the BCS coupling constant and f{E\) is the Fermi function.) As 
usual, the interpretation of the components u\(r) and v \(r) is that they are 
the amplitudes for the quasiparticle being a particle and a hole r espe ct ive!} 



More ove r, a state is superconducting if the solution of equations (2.1), (|2 



and (I2J) is such that A(r) ^ 0. 



In what follows, for a prescribed pairing potential A(r), we develop a 



multicomponent WKB approximation for solving equation (2.1). Asymp- 
totically, as h — > 0, the order of the multicomponent differential equation 
changes abruptly, just as for a single component equation, and therefore, 
we expect the solutions to exhibit the familiar essential singularity in the 
form of the phase e lS ^l h . Thus a reasonable guess at the multicomponent 
generalisation of the WKB wave function would be 

( )-(*«)•***■ (1+otw (gueas) ( "» 

where u\(r) and v\(r) are slowly varying amplitudes. However, to be 
systematic, it is more convenient to start with the general complex spinor: 

u A (r) \_f e^W+sW) 
w A (r) J ~ V e^M- 1 *')) 



(2.5) 



where both S(r) and E(r) are complex quantities to allow for the amplitude 
and phase variation of each of the components. To proceed we expand them 
in powers of h by writing S and X as 

S(r) = ^P- + S 1 (r) + hS 2 (r) + -- - , (2.6) 
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and 



S(r) = 



So(r) 



Ei(r)+ftE 2 (r) 



(2.7) 



In keeping with the original definition of S and E as complex quantities, 
at every order the real (r) and imaginary (i) parts of Sj and Ej must be 
determined. Up to and including the first order in H the spinor in equation 
(p.5[) then contains the following quantities 



e i(E5(r)/ft+SI(r)) e -E|(r) 
e -»(S5(r)/ft+Er(r)) e +S*(r) 



iSS(r)/ft+iST(r)--S , l(r) 



(2.8) 



The question of whether to include Eg (i.e. a fast degree of freedom for 
the spinor components) is a subtle one. If <fi(r) is not expanded in h (and 
it is usual in semiclassical theory to regard potentials such as V(r) and 
A(r) as externally imposed potentials - hence not to expand them) then it 
is found that the BdG equations have no expansion in h if Sg(r) ^ 0. We 
will return to this point again once we have developed our formalism. 
Setting Eq(i") = in ( |2.8| ) we see the wave function does indeed take the 



form (2.4) where however 
varying spinor: 



is understood to be the complex slowly 



u 0) A(r)e ffi i« \ isi{v) e -si{v) 



(2.9) 



where u ,\(r), v ,\(r), and e S ^ T > are amplitudes, and Sq(t), Sl(r) and 
EJ(r) are phases, to be determined by solving the appropriate zeroth or- 



der and first order equations obtained by substituting (2.9) into the BdG 
equations and expanding the result in powers of h. 



2.2. Semiclassical solution of the BdG equations to zeroth 

order 

Noting that e lS °^ r ^ h may be regarded as a unitary operator^ we may 
write, in the position representation, 



-iS (r)/ft- iS (r)/fi = - + 



dSo 
dr 



(2.10) 



Applying this procedure to both sides of equation (2.1) one readily finds 



H(p+^,r) |A(r)|e^ r > / ^(r) 

|A(r)|e-^ r ) -H*(p + ^,r) J { g A (r) 



(2.11) 



1 Since we will not need the r suffix on 5g(r) for what follows we drop it from here 
onwards. 
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Then, the zeroth order approximation in h is obtained by neglecting terms 
containing h (i.e., p = — ifiV, whose action upon u\(r), V\(r) is small). 
Then we observe that 0(r) can be removed from the above equation at 
each r-point by a unitary transformation 

/ „-#(r)/2 q \ 

U4r) e *„ /a • (2-12) 



To zeroth order in h such transformations of equation (2.11) yield 



E x -H§(p ,r) -|A(r)| \ / fi A (r) e -*M/ a \ 

-|A(r)| B A + ^(p ,r) J I £ A (r)e+^ r )/ 2 J U ' 



Here (p , r) is an electron hamiltonian of the form ( |2.2[ ) with the mo- 
mentum operator p replaced by po = ^r, and the same is true for the hole 
Hamiltonian, iJQ(p ,r), except e — > — e. Since the transformed Hamilto- 



nian matrix in equation (2.12) is real, the spinor wave function is also 



constrained. Most generally, via equation (2.£), this implies that 



£r(r)-^=4 (2.14) 

for n an integer. However self-consistency, equation (|j|), requires that the 
A(r) with which we started our calculation, i.e. A(r) = |A(r)|eW, and 
the A(r) constructed from our solutions u\(r) and v\(r) be one and the 
same. Since u\{r)v* x {v) oc e l 0( r )+ m7r 5 we mu st have n = 2m, m an integer, 
for this to be true. Thus we find 

E r {r) _M =mn . (2.15) 

But now, since m introduces the same sign change for both upper and 
lower components of the spinor it can be factored out. Choosing the sign 
of the wave function appropriately m is eliminated and we have determined 

In the context of a general solution of the BdG equations it has been ob- 
served that if the order parameter, A(r), contains vortices, V0(r) acquires 
global curvature p0| i.e 

V0(r)-dr^O, (2.16) 

for paths, c, containing vortices. To clarify the nature of the singularity 
associated with a vortex in three dimensions observe that it is one where 
|A(r)| = along lines. <p(r) is then indeterminate i.e. |A(r)| = defines a 
line of phase singularities |ll| . If we continue ^(r) along any path enclosing 
such a singularity we obtain ( 2.16D with 2ir on the right hand side. A 



unitary transformation of the form is then multivalued, as stressed by 



5 



Anderson p3|. For instance 17$ has two branches and moves from one to 
the other when continued around a vortex. The same is true for our spinor: 
if upon the first branch it takes the values 

then, after a circuit around a vortex, it moves onto the second branch, 
taking the values 

ux(r) 

(A second trip around the vorte x is required to return the spinor to the 
original branch.) Thus, through ( 2.15 ), we have discovered that SJ(r) may 
contain a contribution topological in originQ. This is, however, not a prob- 
lem for our semiclassical theory. With 0(r) removed we must diagonalise 
the resulting matrix subject to the physical condition that the original wave 
function in (2.4) is single- valued. As we shall see presently, this constraint 
will modify the generalised Bohr-Sommerfeld quantisation rule which we 
will derive. 



A non-trival solution to (2.13) requires the vanishing of the determinant. 
This yields two equations: 



£o(Po, 



(2.17) 



where a = ±, and 



£o(Po,r) =Po -vo(r) 




+ m V %(r) + V(r)-e F ) + |A(r)P, 



(2.18) 



where Vo = eA(r)/m. The constancy of i?g(po,r) defines, implicitly, the 
functions Po(r) = 



We see immediately that each of the equations (2.17) is a (stationary) 
Hamilton- J acobi equation for a (fictitous) classical mechanics. Its solution 
(when it exists) is a classical action function Sq(t, I), where Ii,...,I n 
is a set of action integrals and appears in place of A which labelled the 
eigenvalues of the BdG equation. 

Thus viewing each Eq (po,r) as a classical Hamiltonian governing the 
propagation of a quasiparticle excitation the restriction to a constant energy 
shell in phase space (Ef = Eq (po,r)) define the phase space orbits which 
we shall label (po(r, I),r). In general the functions Po(r, I) are many- 
valued functions of r as illustrated in FIG. [|, whose branches shall be 
indexed by j when necessary. Mathematically fl3| finding a complete 

2 It is the singular part of </>(r) which has been called ]lo| a "Berry" phase although 
it is in fact simpler depending only upon the topology not the geometry of the path 
along which it is continued. 
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FIG. 1 Phase space orbits (po(r,I),r) defined implicitly by Ef = 
Eq (po,r). The vector attached to the orbit changes orientation as the 
spinor is carried along the trajectory. Notice that Pq (r, I) is a many- valued 
hmciton of r. 



integral to the Hamilton-Jacobi equation is equivalent to solving Hamilton's 
equations of motion. So we can write down Hamilton's equations for the 
quasiparticle excitation: 



dEZ(p,r)\ /0£ff(p,r) 



dp J ' \ dr 



(2.19) 



where after differentiation we must substitute r = r a (t) and p = Pq (r a (t)). 
Since our Hamiltonian is conservative, time here is simply a parameter 
governing the evolution of a quasiparticle along the orbit. Thus the 2x2 
differential matrix BdG equation ( |2.lD has been reduced, at zeroth order, to 
a pair of Hamiltonian systems for a (fictitous) classical mechanics specified 
by the appropriate Hamiltonian in ( |2.18| ). 



A solution to (2.17) is obtained as 



S^(r,I)= p^"(r,I).dr, (2.20) 

•/ro(to) 

where integration is carried out along a ray of one of the Hamiltonian 



systems (2.19). (For convenience we have chosen the initial condition 
'-'(ro, I) = 0.) This does not complete the zeroth order theory how- 
ever, since our 'particle' has an internal structure which is represented by 
the complex spinor defined at every point along the trajectory. 

So let us now consider the spinor solution of the BdG equation corre- 
sponding to one of these 'classical' orbits. It is useful to view the slowly 
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changing quasiparticle spinor of (2.13) as an axis of pseudo spin quanti- 
sation transported along the trajectory defined by pg(r, I). Clearly, this 
axis can be represented by a vector attached to the orbit (FIG. |l|) whose 
orientation changes as it is carried along the path in phase space. This vec- 
tor represents the internal particle- hole degree of freedom of the excitation. 
To find how the eigenvectors of ( p. 13 ) change as they are carried along the 
trajectory we diagonalise the Hamiltonian matrix, H(r), at every r-point 



(2.21) 



so that in the locally diagonal frame ( [2.13 ) becomes 



E o(Po,r 



£ + (Po,r) 




,,diag 



giving 



E+=E+(p ,r), 



E i = E o (Po,r), 



,diag.+ 



^diag,- 
„diag,- 



■ jdiag 
diag 



(2.22) 
(2.23) 



Rotating back to the common laboratory frame we have two solutions 

u <ti( r ) 
Vi( r ) 



Ej 



e; 



Eo(po,r), 



C/t(r) 
t/t(r) 



Note that while 



is a particle in its local frame, (2.22), it is both 



particle and hole, with amplitudes Ugi(r) and i>oi(r), in the laboratory 
frame. 

If we now substitute Eg into ( p,13| ) we have the ratio 



<j(r) _ £ff-g e ( Po ,r) 
ufotr) -|A(r)| 

The electron Hamiltonian, H{j (p,r) is 

1 



(2.24) 



#o(Po,r) = Po-vq 



Pi 

2m 



V(r)-e F , 



Po-v 



Po-voj 



|A(r) 



(2.25) 



by (2.18), where /? = ±, and the Hamilton-Jacobi equation Ef = Eg(po, r) 
should be used. Using this to eliminate Hg(p, r) from (2.24) the amplitudes 
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are expressed in terms of Eft and the momentum branches, Po(r). It is clear 
from equation (2.25) that each Pq determines a unique choice of (3 (see also 
Appendix |A|) so that there is one pair of amplitudes for each branch of the 
momentum. Then the normalised amplitudes are given by 



u#(r) = 
v oj ( r ) = 



\ 



^i+/ (8f -*'f" 1,11 



Sf-Po-vo 



\ 



1 J(E?-pi- V0 ) 2 -\A(r) 

2 I A £?-P J „-vo 



(2.26) 



(By normalising we have of course redefined e~ s ^ r '^ but since we are yet 
to determine this amplitude we shall not introduce new notation for it.) 
This completes the zeroth order theory. 

To summarise our results so far we have shown that a general semi- 
classical solution for a multicomponent system can be written as a spinor 
multiplying the familiar phase, e lS °/ h (2.4), where however the spinor is in 
general a complex quantity (|2.9D . Using this solution we derived the zeroth 



order matrix equation (2.13) and determined one of the spinor phases, EJ, 
which may have a topological contribution. We then diagonaliscd the ma- 
trix equations reducing the solution of the BdG equations at zeroth order 
to the problem of solving a pair of Hamiltonian systems for our excitation 
specified by the classical Hamiltonians Eq (p , r) . Our excitation has an 
i ntern al structure which, at zeroth order, is represented by the amplitudes 

dH>. 

What is the next step? We expect, to first order in h, to find a transport 
equation whose solution furnishes us with the amplitude ^4(r) = e~ Sl ^ T,I \ 
We also have the phase ^[(r, I) to determine which, like E^(r), is not 
familiar from single-component WKB analysis. Once we have determined 
these quantities our semiclassical wave function will be used to construct 
the rule which quantises our classical dynamics. 



2.3. The transport equation and other first order quantities 

To derive the first order equations we take the expectation value of 
equation ( [2.11 ) and expand the result in powers of h. One can then show 
(see Appendix |b|) that there are two equations to first order in ft: 



0{h) 



1 



Oih) 2Im( u{ v{ 




(2.27) 
(2.28) 
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where P ± — Po(r) ± eA(r). The first of these equations, ( |2.27 ), is the 
transport equation and can be rewritten (see Appendix H) as 



dp 



(2.29) 



P=Po t^ 1 )/ 



It expresses the fact that the product (amplitude) 2 x velocity, at each point, 
must be conserved (there are no sources or sinks). Thus if the local velocity 
of a quasiparticle increases the probability of finding it there decreases, and 
vice versa. Solving equation ( 2.29 ) by the van Vleck JH method gives 



det ^M) 



drdl 



1/2 



(2.30) 



where c is a constant. 

The second equation is more troublesome. It resembles the connection 
enforcing parallel transport in Berry's theory of geometric phases ]l5| , 
p| . However, equation (2.28) cannot be solved to yield a geometric phase. 
Rather we find (see Appendix |b|) 



S[(r) = 



1 



(r) 



V<£(r) 



dt', 



(2.31) 



to 



where 



-r j (r) 



Po J (r) 



+ ((< I >)) 2 -(<'/(r)) 2 ) 



eA(r) 



Equation ( ^.3l| ) is not a line integral. A physical interpretation to S[(r) is 
given in Appendix ||. We should not be surprised that our theory contains 
both a topological phase and another phase which cannot be expressed 
either geometrically or topologically. In a different context, it is known 
Jl7| , pl| p| that the asymptotics solutions of matrix differential equations 
can contain such terms. 

Bringing all our results together, the multicomponent semiclassical wave 
functions for the Bogoliubov-de Gennes equations with action Sq' 3 (r, I) 
take the form 



u?(r) 
wf(r) 



det 



drdl 



1/2. 



u$(r)e+**M/ 2 
»^(r)e-**W/ 2 



exp^fi-^Cr.IJ+t^'ir,!)) , (2.32) 



where is a constant. 
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2.4. Torus quantisation 

To construct a quantisation rule we must consider the behaviour of 
( ||§ ) along rays of the Hamiltonian system (|2.19|). Thus we must first 



discuss the properties of the underlying classical dynamics. 

It is well known p(J that for a classical mechanics with N degrees of 
freedom to be integrable there must be N constants of the motion. For 
such a system the dynamics takes place on an iV-dimensional (Lagrangian) 
submanifold in phase space, which has the topology of an TV-torus [^3). To 
perform a canonical transformation from old coordinates (p, r) to new co- 
ordinates (I, <p), where ip are angle coordinates on the torus, one constructs 
the action function, S^r,!), which is a solution of the Hamilton- Jacobi 
equation for the given classical mechanics. Then the transformation is 
specified by 

dS(r,I) dS(r,I) 

The action variables, I, which are constants of the motion, are given by 

1 



pdr, (2.33) 

2vr J Tl 

where Ti is the Zth irreducible loop on the iV-torus. The energy can then 
be expressed solely in terms of these constants, E — E(T), and Hamiltons 
equations of motion can then be integrated explicitly. 

It is also well known j^, |7) that the semiclassical approximation to 
the spectrum of the Schrodinger equation with a Hamiltonian whose clas- 
sical dynamics is integrable, proceeds via the semiclassical wave function 



(whose multicomponent generalisation we have given as (2.32)). The single- 
valuedness of this wave function when followed around each of the irre- 
ducible loops upon the torus yields the N quantisation conditions 

p • dr = 2nli = 2irh (ni + ^j-) , (2.34) 

where ni and m; are integers, mi is the Maslov index ||, and accounts 
for a change of n/2 in the phase of the wave function each time a caustic 
is crossed. The Maslov index of a closed curve is defined to be the number 
of times dr/dp changes sign from negative to positive, minus the number 
of times dr/dp changes sign from positive to negative. It is a topological 
invariant of a Lagrangian torus. 

The quantisation conditions ( [2.34 ) are the so called Einstein-Brillouin- 



Keller (EBK) quantisation rules, the generalisation of the Bohr- Sommerf eld 
quantisation rules of 'old' quantum mechanics to encompass non-separable 
problems. 

What we require is the multicomponent generalisation of the EBK rules. 
Since the multicomponent theory was reduced to two Hamiltonian systems 
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one might think that the EBK rules apply to the tori for each Hamiltonian 
system. We start by writing the most general semiclassical wave function, 
which due to the multivaluedness of Pg(r) is given by the superposition 
principle of quantum mechanics as 



uf(r) 
itf(r) 



E 



det ^M) 



x exp 



drdl 



1/2. 



X 0,I 



(r)e 



+#(r)/2 



(r,I) 



-iS"' 3 (r,I) + *nwr/2) , (2.35) 



where j runs over all the branches, Sq , which together make up the 
torus. Firstly consider a situation where S^' 3 = 0. Such a situation occurs 
when = 0. Then it is clear that single- valuedness of this wave function 
does yields a quantisation condition of the form ( 2.34 ). Now generalise 
to the case where V<^(r) has global curvature. Then the contribution to 
the change of phase around a closed path arising from the topologically 
non-trivial <fi(r) is wm^ where 



= — 4 V.Xrl • </r. 



(2.36) 



is an integer. However the phase S"' 3 (r) poses a greater problem since 
as discussed it is not locally path independent so that it prevents us from 



constructing an analog of (2.34) in the most general circumstances which 
we will be interested in. This problem has been discussed in the context 
of multicomponent wave equations by Littlejohn and Flynn 17 . They 



proposed using Hamiltonians which include first order terms in h. For these 
new Hamiltonians the effect of S"' 3 is already included in the 'zeroth order' 
theory. Their results were limited to systems with no global degeneracies 
and, as discussed by Emmrich and Weinstein Jl9|], integr ability of the 
underlying classical dynamics does not garantee an extension of the EBK 
rules for degenerate systems. Recently Bolte and Keppeler (2l| presented 
a semiclassical theory for the Dirac equation based on semiclassical trace 
formulae, rather than multicomponent WKB, at least in part to avoid these 
difficulties (and also to handle non-integrable dynamics). In our case, at the 
present time, this will not be necessary. In the next section we will follow 
a similar proceedure to Littlejohn and Flynn which however will differ in 
that the Hamiltonians we construct will naturally be seen to contain terms 
up to h 2 . (This takes our work outside of the considerations of the above 
citations |ll| |l7||.) Later in this paper, when we study the single vortex, 
the inclusion of h 2 terms in the Hamiltonian turns out to be essential in 
order to obtain the known (exact) wave function at the origin. (Note 
that there is no reason to object to ^-dependent terms in the 'classical' 
Hamiltonians since it is an effective (fictitous) classical mechanics that we 
are considering.) 
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3. AN EFFECTIVE SEMICLASSICAL THEORY 



3.1. The zeroth order theory 

We start by writing the spinor in terms of amplitudes and phases 



Mo ,A(r)e+ lS M 
« ,A(r)e- 4S ( r ) 



e iS(r) e -S;(r) 



(3.37) 



but this time we don't expand S and £ in h. The amplitudes uo,a(f), fo,.\(r) 
and e~ Sl ^ are again determined by the zeroth and first order equations 
of the theory but these equations will be different from those obtained by 



asymptotics in H. Substituting (3.37) into the BdG equations and proceed- 



ing as before we obtain, instead of (2.11), the equations 



-|A(r)|e-^W 



-|A(r)|e^« 
-tf*(p + p(r;7i),r) 



x e 



u a .x(r)e+^ 
v ,x(r)e- l ^ r) 



0, 



where p(r; h) = hS7 S is the S-dependent momentum^]. We are not yet ready 
to make any simplifying approximations since the spinor still contains a 
phase. We use U^{r) (2.12) to remove <p(r) at every r-point and obtain 



E x - H(p + p(r; H) 
-|A(r)| 



fV0,r) -|A(r)| 
Ex + H*{p + p(r;h) 

u ,x(r)e+^~^/2 \ _ sUt) 
v Q .x{r)e-^+^/ 2 ' 



2 v 



r) 

0. (3.38) 



For an appropriate zeroth order theory the spinor written in this form is 
indeed real as we now demonstrate. 

To simplify the equations ( [3.38] ) we must drop the differential operators 
whose action upon uo t \(r) and vo t \(r) is small. Previously this was achieved 
by discarding all terms containing h (p — — iSV). However consistency 
would then require us to replace p(r; K) by po(r) which we do not want. 
Rather than use H as our ordering parameter we instead use p itself. Thus 
our zeroth order theory is obtained by dropping p to give the analog of 
O i.e., 



Ex 



H e (p,r) 
|A(r)| 



Ex 



|A(r)| 
H h (p,r) 



U A (r) e +iS(r)-i0(r)/2 

v ,x(r)e~ tnr)+l4,{T}/2 



= 0, 



where now p = p(r;h), H e (p,r) is an electron Hamiltonian of the form 
( [2.2| ) with the momentum operator p replaced by p(r;?j), and the vec- 
tor potential A(r) replaced by an effective vector potential A cff (r) = 

3 If we wish to recover our previous results we can expand: p(r; K) = VSo + ^VS^ + 
■ ' ■ = Po(r) + ^Pl(r) + ■ ■ ■ , and proceed with asymptotics in the small parameter h. 
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5jV0(r) + A(r). The hole Hamiltonian, H h (p, r), is the same as H e (p, r) 
except e — > — e. 

Since the resulting matrix Hamiltonian is real we have 

E(r) = (3.39) 

and again, as a consequence of ( [2.16 ), we see that £ can contain a topo- 
logical contribution. 

Our new theory yields the Hamilton- Jacobi equations 

£i = £ Q (p,r), (3.40) 

where two new /i-dependent Hamiltonians, E a (p, r), have replaced the 
Hamiltonians Eq (po,r) appearing in our previous theory. 
Explicitly the new Hamiltonians are 




E a (p,r) = p • v s (r) + a y + ~mv 2 (r) + V (v) - e F j + |A(r)P, 

(3.41) 

where a = ± and tov s = |V0 + eA(r). The constancy of _E Q (p,r) defines 
implicitly the functions p a (r; h) = , and the solutions of the Hamilton- 
Jacobi equations, when they exist, have the form 

S a {v,l) = \ p a (r,I;H)-dr, (3.42) 

n Jr (t ) 

where the integrals are taken along trajectories of the Hamiltonian systems 
fdE*(p,i)\ fdE a (p,r)\ 



\ dp ) \ dr 



The //.-dependent Hamiltonians (3.41) together with the Hamiltonian sys- 
tems they define, and the ^-dependent action (3.42) are a central result of 
this paper, ^"(r, I) is clearly a locally path independent quantity since 
p a (r, I; fi) is a solution of Hamilton's equations and thus lies in a La- 
grangian torus in phase space. If we can prove that no other phases appear 
in our theory the action can immediately be used to construct a quantisa- 
tion condition. 

You will notice that the Hamiltonians (3.41) include a term in (r), 
that is a term of order h 2 as alluded to in the preceeding section. 

In concluding the zeroth order theory we see that the normalised am- 
plitudes, representing the internal structure of our 'particle' as it is trans- 
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ported along the trajectory, become 



0,1 V-) - \ 2 I i_r ^ E°>-pl- Vs 



1 A _ ^V(-E°-p J -v 3 ) 2 -|A(r)| 



<i (r) = y § ^ - /? v ^ xv/.v7 — J ■ ( 3 - 44 ) 

Thus the amplitudes too have become ^-dependent. 

3.2. First order theory 

To derive the new first order equation we take the expectation value of 
equation ( 3.38) ) . By defining 

uo,A(r)e+ ffi W-*W/2 \ „ s;(r) 



G= { voT(r)e-^r) + i< t >(r) / 2 ) e~™, (3.45) 

and introducing D(h) for the ^-dependent matrix differential operator: 

- / ^-^(p + p^^ + fV^r) -|A(r)| \ 

[ ) \ -|A(r)| £ A + J ff*( P + P (r;ft)-fV0,r) )' 

we can write this expectation as 

G*D{h)G = 0. 

Expanding this equation upto and including first order in p we find 

= GlD Q (h)G + G\D Q (h)G + GjA, + Gj£>i(fi)G , (3.46) 

where -Do (ft) is the zeroth order ^.-dependent Hamiltonian matrix 

nm_/fi-ffV) -|A(r)| \ , , 

Do(S) -i -|A(r)| Ei + fl*(p,r) )' (3 ' 47) 



Go is given by equation (3.45) together with the zeroth order condition 
( p9| ) i.e. 

Go=f ^'tn e" s H (3.48) 



«o,i(r) 

which does not contain any phases, Gi allows for corrections beyond those 
considered in Go (in analogy with the semiclassical theory (Appendix |b|)), 
and D\{h) is 



Di 



-2k(iV-P + + P + -|v) \ 
o +^(fv-p-+p--fv) J' 
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where P ± (h) = p(r; h) ± mv s (r). Our analysis now parallels our previous 
semiclassical discussion (Appendix |b|) : the first term in equation (3.46) is 
the zeroth order equation whilst the second and third terms can both be 
shown to depend upon Do(h)Go = 0, so that their vanishing yields no new 
information. We are left with the first order equation 

<D(p) = GjDi(S)Go. 
Defining the vector matrix M by 



M = 



-™ 

2m 



2m 



this can be rewritten as 



= V 



{g o MG(T 



2i Im 



GjM • 



VG 



(3.49) 



(3.50) 



The first term is purely real, whilst the second is purely imaginary. In our 
previous semiclassical theory, where we had the spinor Fa in plac e of Go, 
the first term led to the transport equation (see ( [2.27 ) and (|2.2S| )) whilst 
the s econd led to an equation for the first order phase S{(r), (see ( 2.28 ) 
and (2.31)). This was because F was complex. In the present effective 
semiclassical theory Gq is real ( [3.48 ) . Consequently the second term in the 



above is trivially zero, 
can be written as 



We are left with a new transport equation which 



\ dp 



= 0, 



(3.51) 



in terms of our new ^-dependent Hamiltonians ( 3.41 ). Its solution yields 
the new /i-dependent determinant 



-Sl(r,D = c 



det 



d 2 S Q J(r,I; ft) 



drdl 



1/2 



(3.52) 



Actually one must be more careful than we have been here in deriving 
this effective theory. Whilst one would very much like there to be no first 
order phase corrections (S[, S^) we simply cannot demand that this is the 
case. As soon as we approximate to obtain a zeroth order theory S and S 
are only determined to this order and we must allow for the possibility of 
corrections, ^[(r), S^(r). If we do this G, equation (3.45), is replaced by 



r<' — [ "O.IV 1 ^ I p +iS{ r) -S\(t) /o r,o\ 

_ Uoi(r)^ lS(rhlE;(r)+ ^ (r)/2 J ' ^ > 
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and Gq by 



a' — ( u o,i(r)e+* s i» ^ r < 



(3.54) 



Then Im 



GqM • VG 



is no longer trivially zero and yields 



where 



— J ( r ) = - 

e 



Sl(r)=-- /" j^(r).VEl(r)< 

+ ((<;/(r)) 2 -(<; i '-(r))^ " A " j!fr) 



(3.55) 



'•(r;R) 



But there is a crucial difference between (3.55) and (2.31), namely it de- 
pends upon VSJ rather than V</>. What can we say about VEJ? Appealing 
to self-consistency we have 

UA (rK(r) OC e *2£(r)+^(r) = e #(r)+,2EKr)^ 

Thus = rmr, to an integer, so VE^ = 0, and hence S[ = 0. Furthermore 
using the same arguments as followed equation (2.15), the need for Y,\ 
can be eli mina ted. Since S[ = 0, and = o ur form for the spinor, 
equation ( [3.37 ), and all that follows up to ( 3.51 ) is indeed correct. We 
have succeeded in constructing a theory where there are no first order 
phase corrections and have thus removed the obstacle to the derivation of 
a generalised Bohr-Sommcrfeld or EBK quantisation rule. 

The general wave function for our effective semiclassical theory now 
takes the form 



uf(r) 
i>f(r) 



det 



d 2 S a 'i(r,I;h) 



drdl 



1/2/ 



u o,i,Jl( r ) e " 



x exp (ih^S^ir,!; h) + imir/2) , (3.56) 



and the single- valuedness of this after returning from circuits around each 
of the irreducible loops, Ti on the 3-torus yields the general quantisation 
conditions 



<j> p Q (r;7i)- dr = 2nh 




where as before 



(3.57) 



(3.58) 
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takes integer values. For a given T; the line integral in equation (3.57) 



is fixed so we must have = n, + mf, i.e. we need only one of nr 



Choosing nj~ = n\ equation (3.57) becomes 



p a (r;/i).dr = 27rftfn I + ^-^j ) (3.59) 



for ni, mi, and mf integers. This is our generalisation of the EBK quan- 
tisation condition to apply to the superconducting case. It includes two 
topological integers. The first is the familiar Maslov index, whilst the sec- 
ond arises from the vortex singularities, which, if present, shift the quantum 
numbers by half- integers. The action integral itself, as discussed, is defined 
upon an ft-dependent Lagrangian submanifold in phase space. The tra- 
jectories of the Hamiltonian system which wind around this manifold are 



specified by the ^-dependent Hamiltonians, E a (p, r), equations (3.41). The 
question then arises as to the interpretation of a theory depending upon 
these Hamiltonians. In particular how should the appearance of hV<fi(r)/2 
be understood? 

3.3. Interpretation: A Semiclassical theory in the presence of 
lines of phase singularities 

Consider the BdG equations written in the (exact) form 

( ^A-^( P + fV0(r)+ e A(r)) 2 -nr) + e F -|A(r)| \ 

^ -|A(r)| E X + ± (p - fV-Kr) - eA(r)) 2 + V(r) - e F J 

UA ( r ) e -#(r)/2 



x ( : A ^U( r )/ 2 ) = 0. (3-60) 



The term fiV0(r)/2 enters into the particle and hole Hamiltonians as an 
effective vector potential 

A PS (r) = A V0(r). (3.61) 

ze 

There is a flux associated with this vector potential which we can find by 



integrating (3.61) along a path enclosing a vortex singluarity: 

j( V</>(r) • dr = A j( a ps (r) • dr, 

2tt = A$ Q . (3.62) 
n 

The left hand side follows from the single- valuedness of A(r). On the right 
hand side we have introduced $o for the flux. It takes the value 

$o = — , 3.63 

e 
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i.e it is equal to one superconducting flux quantum. Where is this flux? 
We can contract the curve, c, until only the line of phase singularities is left 
inside. Equation ( 3.62| ) remains true. Thus each line of phase singluarities, 



defined by |A(r)| = 0, carries a flux, $ , along it. A similar object was 
studied by Dirac ^2j, and is refered to as a Dirac string. However our 
string has a physical reality, lying along the node of the order parameter 
unlike a Dirac string which need not lie along a node. Hence we will stick 
with the 'line of phase singularities' terminology. A PS (r) is then the vector 
potential associated with such a line of singularities. We can then interpret 



equation ( 3.6C ) as describing superconducting quasiparticles in the presence 
of lines of phase singularities of the pairing potential, each carrying a flux, 
$o- Replacing p by p and diagonalising we obtain the Hamilton- Jacobi 



equations E\ = E a (p,r), with E a (p,r) given by equation ( p.41[ ). Thus we 
have constructed a (fictitous) classical mechanics describing quasiparticle 
excitations in superconductors propagating in the presence of lines of phase 
singularities carrying flux $o = h/2e. 

When the classical mechanics is integrable (a single s-wave vortex is such 



an example) then the semiclassical wave function takes the form ( 3.5q ) and 
we can apply the quantisation conditions ( |3.59| ) to obtain the semiclassical 
excitation spectrum. If however the classical dynamics is non-integrable 
then the solution cannot have the form ( |3.56| ). We return to this point at 
the end of the paper. 

In the interest of clarity concerning the above discussion we would like 
to make the following remark. If for a moment we consider a single vortex 
whose axis lies along the z-axis, then by symmetry 

A(r) = |A(r)|e ie , 

in polar coordinates. The vector potential A PS (r) is then explicitly 

A PS (r) = — 6. (3.64) 
2irr 

Now the vector potential associated with an Aharanov-Bohm (AB) flux 
tube is 

<s AB 

AAB(r)= 2^' (3 ' 65) 
for r > r c , the solenoid radius. If we consider the limit r c — > as repre- 



scnti ng an idealised AB flux tube, or AB flux line then equations (3.64) 



and ( 3.65 ) become identical for the choice <I> AB = $o- For this reason a 
vortex has been described in the literature [jioj] as an effective magnetic 
Aharonov-Bohm half-flux. However we would like to point out here that 
there is a crucial difference between an AB flux line and a superconducting 
vortex, namely that the single-valuedness of the AB wave function does 
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not quantise the flux, <E> AB , which remains classical ^3|, whilst single- 
valuedness of the order parameter, A(r), does quantise the flux, <I>o- For 
this reason we prefer the line of phase singularities terminology. 

We have one last point to make. In our original discussion of the general 
form for the wave function we remarked (see the paragraph following equa- 
tion (2.8)) that if <fi(r) is not expanded in h, which we claim to be the correct 
procedure, then there is no fast component to the spinor i.e., E(j(r) = 0. 
Conversely had we insisted upon expanding (f>(r) = ftr 1 ^^) + ■ ■ ■ the ze- 
roth order theory would have included the term fiV(f>\r l= o — V(fo- In the 
light of our effective semiclassical theory we now see that even though </>(r) 
is not to be expanded, i.e., ftV</>(r) ~ H, the most useful form for the theory 
has fi.V(/> appearing in the "zeroth order" theory Hamiltonians. 

This concludes our formalism for seeking semiclassical solutions to su- 
perconducting problems described by the BdG equations (^j]) . The general 
procedure is then as follows: (i) Choose A(r) and A(r) appropriate to the 
situation of interest. One can then immediately write down the effective 
Hamiltonians in explicit form, (ii) Investigate the orbits of the Hamiltonian 
system and identify those of interest, (iii) Compute the relevant Maslov 
indicies and any topological indicies arising from vortex singluarities. (iv) 
Deploy the generalised EBK quantisation condition ( [3.59 ) to obtain the 
semiclassical spectrum. (Steps (iii) and (iv) can be interchanged if analytic 
expressions are being sort.) Depending upon the situation it may also be 
important to: (v) investigate tunnelling between classical orbits and its 
effect upon the semiclassical spectrum, for example the lifting of degen- 
eracies to yield 'avoided crossings', (vi) One may also deploy the wave 
function to study properties other than the energy spectrum (for example, 
quasiparticle current flow). 

In the next two sections we will apply this theory to two well known 
problems. The first we will consider is a superconductor-normal metal 
-superconductor (SNS) junction. For this system we will focus upon un- 
derstanding the various branches of Py(y) which make up a trajectory in 
a superconducting system, and upon the type of quasiparticles they repre- 
sent. We will then quantise the orbit and obtain the spectrum. Further- 
more we will compute the Maslov index, though via the more 'traditional' 
route of asymptotic analysis, and show that it is what one would expect 
from following Maslov's topological prescription. Since for the SNS junc- 
tion we take hW(j> — the two theories presented in sections || and || are 
identical. By contrast the second problem we apply our theory to, the 
single vortex, has HS/(j) ^ 0. The significance of the topological phase and 
ft-dependent terms in the theory will then become apparent. 



20 



4. THE SUPERCONDUCTOR-NORMAL 
METAL-SUPERCONDUCTOR JUNCTION 



A superconductor-normal metal-superconductor junction (SNS) con- 
sists of a taking a superconducting wire and replacing a segment with a 
piece of normal metal. The inclusion of the normal layer makes the order 
parameter, A(r), inhomogeneous along the wire. By assuming the wire is 
sufficiently thick to neglect finite size effects the profile of A(r) will only 
vary along the length of the wire. We choose the y-axis to lie along this 
direction and take y = to be at the middle of the normal layer. We shall 
allow |A(y)| to have a smoothly varying profile at each interface. This 
contrasts with the more usual approach in the literature |2[| of taking 
the profile of |A| to be a step function at each interface. The 'width' of 
the normal layer, as far as an excitation is concerned, is decided by the 
turning points of the classical trajectory and is thus energy dependent. 
These turning points will be designated by y-(E) and y+{E) (y_ < y+), 
see FIG. 0. The phases of the order parameter are taken to be constant and 




FIG. 2 Illustration of the profile of |A(y)| with classical turning points, 
y±(E), indicated. 



both equal to <j>, i.e there is no phase gradient across the junction. Then, 
taking the vector potential A(r) to be zero, the classical Hamiltonians, 



■E a (p,r), equations (3.41), are found to be 




B a (p, r) = ad + + £- - e F ) + \ A(y)\*, (4.66) 



where a = ± distinguishes the two Hamiltonians. However in keeping with 
the definition of E\ as an excitation energy (E\ > 0) we can discard the 
a = — 1 Hamiltonian, and consider only E + (p,r). 



21 



4.1. The y behaviour of the excitation orbits 

Setting E + (p,r) — E, equation (4.66) can be inverted giving p^(y) 



(/? = ±) as 

pP(y) = yfe -pl-pl+ (32m^Ei - | A(^. (4.67) 

Thus E + (p,r) = E defines four momentum branches (including — jc(j/)) 
as a function of position, from which the orbit of the excitation is to be 
constructed. 

Now Py represents a quasiparticle, and p~ a quasihole (see Appendix 
^). This is clear from considering (4.67) deep inside the normal layer where 
[A(y)|-0: 



Km p±(y) = s Jp%±2mE-pl-pl (4.68) 

When E is non-zero p+ lies outside the Fermi sea - a particle- like excitation 
- whilst p~ lies inside, so describes a hole-like excitation. 

Throughout the normal region, for which |A(y)| = 0, the momentum 
branches Py (y) are independent of y, i.e are straight lines in the p y — y 
phase plane 

i ■ 1 1 1 ' i , 



Let us write Py(y) in (4.67) as 



tf(y) = \ p 2 f ~pI-pI + Hv), ( 4 - 69 ) 



where e(y) = 2m\/E 2 — |A (?y)P . Approaching the interface |A(y)| — » E, 
so e(y) — > + . Then from ( 4.69 ) p^ decreases and p^ increases (FIG. ||). 



The length scale over which e(y) changes is given by the coherence length, 
£. At the interfaces e(y±) = and Py(y) are 



lim p±(y) = Jp 2 f -pI~p 2 z . (4.70) 

|A|— >_fc v 

Thus at the points y±, the particle and hole momenta are equal. Notice 
however that p y z (y±) 7^ 0.^ This behaviour at the coelesence of the mo- 
mentum branches contrasts with typical classical orbits where the turning 
points are characterised by p y {yi) = 0. 

To understand what happens when the two momenta coalesce we con- 
sider the velocity associated with an excitation moving along each of the 



4 From equation (4.70), p y (y±) ^ unless p 2 F = p^,+P^, i.e. unless all the kinetic 
energy is parallel to the interface - a situation which we shall not consider here. However, 
see for example Sipr and Gyorffy for details of how the physics changes radically 

in this limit. 
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FIG. 3 The form of Py{y) as |A(y)| — * E. The p^ are constant and 
positive deep in the normal layer y_ -C y -C y+. For |A(y)| — > E 1 (y — » j/±) 
p+ decreases whilst p~ increases. 



trajectories p^. F or th is we use Hamilton's equations ( 2.19 ), together with 



the Hamiltonian ( 4.66 ) for this problem. We have 

(dE(jp, r y 



r = 




p 

2m 



eF + |A(l/)| s 



?7) 



all evaluated for p 
velocity is 



p ± (r) and r = r ± (t). By using equation (4.67) the 



= ± 



y/E*-\A(y)\* (p x ,P±(y),Pz) 
E m 



(4.71) 



In the limit |A(y)| -> 0, v± -> ± (p x ,pf(y),p z ). Let p± > as in FIG. | 
then for a particle the velocity v + is directed to the right, whilst the hole 
velocity is directed to the left. This situation is reversed if the negative 
branches of p^ are considered. FIG. || shows all the branches of p y (y) 
together with arrows indicating the velocity of the excitation. Consider 
a particle travelling along p+ > approaching y + . The velocity v + > 
decreases until it is zero at y + . At this point p+(y+) = p~{y + ) and the 
particle converts to a hole. It then moves away from the interface with 
v~ < 0. Thus even though Py(y) ^ 0, y — y± is a turning point of the 
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FIG. 4 The branches of p y (y) with the corresponding velocity directions 
indicated. 



orbit. From (4.71) we see the quantity responsible for the change in sign 
of the velocity is the scalar function ±y/E 2 — \A(y)\ 2 /E which multiplies 
all three components of the momentum p. Thus at reflection all velocity 
components reverse their direction despite both p x and p z being constants 
of the motion.^ Such a process is called Andreev reflection ^ , and we have 
demonstrated that our 'classical' Hamiltonians describe this remarkable 
property of quasiparticle excitations in a superconductor. 

Once the hole has travelled across the normal layer the reverse process 
happens, the hole turns into a particle, and the excitation completes a 
periodic orbit. Thus we have bound excitations in the normal layer whose 
orbit consists of both quasiparticle and quasihole segments. We call y± 
(defined implicitly by e(y±) = 0) Andreev turning points to distinguish 
them from normal turning points, y%, for which p^j (jji) — 0. FIG. ^conveys 
this remarkable feature of the reflection in real space. 

Now that we have identified the classical orbits we can turn our atten- 
tion to the semiclassical spectrum. 

4.2. Calculation of semiclassical spectrum 

We deploy our generalised quantisation rule 

f rn rri W 

V a {r;h)-dr = 2-Kh[n+ 



5 It is clear that p x and p z are constants of the motion because both x and z are 
cyclic coordinates so that p' x = and p' z = 0. 
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FIG. 5 Illustration of the unique nature of Andreev scattering - retrore- 
flection. 



Since p x andp z are constant along the orbit, which retraces itself, § p x dx — 
§ p z dz = 0. Furthermore to^ = since there are no vortices. Our quanti- 
sation condition becomes 

£pP(y;h)dy = 2nh(n+^, (4.72) 

with = ± along the appropriate segments of the t rajectory, T, drawn in 
FIG. f|. If we calculate the right hand side of ( 4.72| ), which depends upon 



E through pP, we can find the semiclassical spectrum of bound Andreev 



excitations in the SNS system. This we now do. 
For our example we have 

P P y (V) dy = f + P + (y) dy - f + p~ (y) dy. (4.73) 

Jy— Jy- 

At this stage we approximate the momentum branches Py (y) by replacing 
|A(j/)l by a step |A|. There is no need to do this[] other than to facilitate a 
comparison of our results with Andreev |3J . With this approximation the 



iix O, A 

fingfA 



In appendix tj, when the Maslov ir 



smoothly varying^A(j/)|, see equation (C.142) 



Hex rp .i is calculated, it is done so for a general 
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momentum branches p^ in the normal layer are given by 

Py = \]P F -P±± 2mE > 

2mE x 1/2 



V P f ~ P± 

where p\ — p 2 x + pi- The quantity 2mE/(p 2 F — p\) <C 1 for almost all 
values of p± except p± — > pp. The latter limit implies that the excitation 
hits the interface at glancing incidence which we do not consider here. Let 
us therefore expand Py in powers of E. Upto and including first order in 
E the p^ become 

+ ■ „, , 1 2mE 

p±=p F cos6 ±- — -. 4.74 

y 2p F \cos0\ 

Here we have used p± — p F s'md (valid to this order). (The modulus has 



been inserted to remind us that — 7r/2 < 9 < n/2.) Returning to (4.73) 



and using (4.74) we have 



y+ 



2E 



dy (Pl(y)-Py(y)) = " a . L, (4.75) 

y " u_f|cos£^| 



where L = y + — y_ is the width of the normal layer which is the same 



for all bound excitations in the step |A| approximation. Quantising (4.75) 



directly using (4.72) we find ||: 

E n (e)=nHv F (n+^)^. (4.76) 

This is the well known spectrum of an excitation trapped in the normal 
layer of an SNS junction (usually obtained by wave function matching). 
Since the orbits we have qua ntised are topologically circles the Maslov 
index is m = 2 (see section |2.4| ). We will however ca lcula te it explicitly be- 
low. A number of text books discuss the spectrum, ( [4.76 ), (see for example 



Abrikosov |27]]) so we will not dwell on it further. We do, however, want 
to draw attention to the following points. Our use of Hamiltons equations, 
together with a smoothly varying |A(y)|, gave a simple and clear picture 
of the 'particle' to 'hole' conversion at the interface. It is an attractive fea- 
ture of this theory that our classical Hamiltonian system describes Andreev 
retroreflection, rather than having it arise from wave function matching. 
The generalised quantisation rule correctly reproduces the Andreev spec- 
trum. In his book |27| Abrikosov invokes Bohr's quantisation rule to obtain 
this spectrum. However the justification for using such a quantisation rule 
is non-trivial as we have demonstrated. 

We now conclude our discussion of the semiclassical theory in the con- 
text of the SNS junction by utilizing the wave function to determine the 
Maslov index. 
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4.3. The semiclassical wave functions for excitations in SNS 

junctions 

The wave function for this example consists of plane waves in the x 
and z directions (since p x and p z are constants of the motion) together 
with the one dimensional form of our semiclassical wave function. (Since 
V0 = either of the semiclassical wave functions, ( [2.32j ) or ( |3.56| ) in their 
one dimensional form can be used because they are identical in this case.) 
Thus we have 



0,3 



0,3 



d 2 S^{yJ) 



dydl 



x l 2 ( „,0,j 



.,0,3 
V 0,I 



i(f>/2 



x e ipxX ' h+ip ^ n exp (ih^S^iy,!) + irrm/2) , (4.77) 



where now (3 labels the two distinct momentum branches Py and p y whilst 
j specifies the sign of those branches. It is useful to separate the four 



branches Py , p y 



-Pi 



-p y , in this way because both the spinor amplitudes 



u a'i(y) ana - v o'i(y), an d the amplitude factor \d 2 S^^ / dydl] 1 / 2 will turn 
out to be insensitive to the sign of the momentum. The amplitude factor 
can be rewritten as 



d 2 S^(y,I) 



dydl 



1/2 



1/2 



dE 



dp y 



-1/2 



(4.78) 



where — dE/dl is a constant which we will absorb into the ampli- 
tudes A 0jj . We then have S^(y,I), u^(y), v^'j{y) and \dEjdpy\- 1 ' 2 to 
investigate. 



4-3.1. Spinor amplitudes 



The spinor amplitudes are given by ( 3.44 ) with v s (r) = and |A(r)| 
\A(y)\ i.e 



u oj(y) 



\ 



1 / /3y/E* - \A(y)f 

2 I E 



(4.79) 



and 



v oj(y) 



\ 



1 ^E 2 -\A(y)\< 

2 I E 



(4.80) 



Equations ( 4. 79] ) and (4.80) show that the semiclassical expressions for the 
spinor amplitudes reproduce the familiar form for the SNS junction problem 
(see for example f28|). Briefly consider the two limiting cases |A(y)| — > 0, 
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and |A(y)| — > E. In the first instance, deep inside the normal layer, the 
spinor amplitudes corresponding to Py and p~ are easily seen to be 

l™ ( ^/fl ) = ( l), (particle), 



|A(y)|-M3 V V 0,l(v) J V 

Km ( ) = f J V (hole), 

|A(»)|-0 V «0,l(») / V 1 / 

so that the spinor amplitude for the particle-like excitation, corresponding 
to the momentum branch pt(y) is (J) , and for Py(y), it is (J) . Notice 
then that particle- and hole-like excitations decouple in the normal layer, 
as is to be expected. An excitation is either particle or hole, not a mixture. 
In the second limit, |A(j/)| — > J5, we have 

«m ( ) . f 7 V 

Hm f ) = f ? y 

Thus at the turning point an excitation is an equal mixture of particle and 
hole. If we assign an effective charge e*(y) — e (v,q j(y) — Vq I (y)) to a given 
spinor then a particle- like excitation has e*(y) > 0, a hole- like excitation 
has e*(y) < and at the turning point e*(y) = 0. (Recall that charge is 
not a good quantum number in a superconductor.) As a final comment, 
notice that for |A(y)| > E, Uq T (y) and v$ j(y) become complex conjugates. 
This coincides with the momenta Py(y) becoming complex. 

4-3.2. The velocity dependent amplitude 



The wave function (4.77) contains not only the position dependent par- 
ticle and hole amplitudes j(y) and Vq j(y), but also the overall amplitude 
(dE/dpy 1/2 . We have 

f dEjy^ y 1 ' 2 £ 1/2 mV 2 



Here the y/Py{y), familiar from semiclassical theory for one component 
systems, is supplimented by another term. For the SNS junction we again 
consider the two limits |A(y)| — > and |A(y)| — » E. In the first instance 
we have 



lim 



fdE(p,y)Y 1/2 m 1 / 2 



|A(v)HoV dp y ) Vy=vl{v) 
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so that a particle-like excitation has an amplitude inversely proportional 
to the square root of the particle momentum, and the hole-like excitation 
similarly depends upon the hole momentum. In the second limit as the 
excitation approaches the interface (dE/dp) -1 ^ 2 diverges, despite Py{y) ^ 
0, due to the (E 2 -\A( y )\ 2 )^ 1/i dependence. This divergence signals the 
breakdown of the semiclassical approximation for the wave function when 
we are too close to the turning points (i.e., caustics) at which the velocity 
goes to zero. 

4.3.3. The action S$' j 



The action appearing in (4.77) is given by the integral of the appropriate 
momentum branch. We have 



S^(v)=j F dy'p%(y'), (4.82) 

Jyo 

where j = ± gives the sign of the momenta, and (3 as usual distinguishes 
the branches Py and p~ . The lower limit yo is an arbitrary constant called 
the phase reference point. 

Although the amplitudes Uq r (y), Wq r (y), and (dE/dp)^ 1 ^ 2 acquire 



their normal state forms as |A(y)| — > 0, the action (4.82), being an in- 
tegrated quantity, retains a 'memory' of the interface region encountered 
by the quasiparticle. 

If we consider the form of Sq' 3 (y) inside the superconductor (|A(y)| > 
E) then the momenta, pf^(y), become complex and the action correspond- 
ingly has both a real and imaginary part: 

So - 3 (y) = i f d y' pW) + Hi/?) f H M), (4-83) 

where p r y and p y are the real and imaginary parts of the momenta. 

4.4. Derivation of the quantisation condition and 
determination of the Maslov index 

We now construct a specific solution. In all that follows we will con- 
centrate upon the y-dependence of the wave functions (i.e., we will drop 
the plane wave factors for the x and z directions since they play no role 



in this derivation). As a consequence of (4.82) the wave function in the 
superconducting region is comprised of only those solutions which satisfy 
the physical boundary conditions, decaying as y — > ±00. Let us label the 
three distinct regions of the SNS junction A, B, C FIG. | where B is the 
'classical' region t/_ < y < y + . Then in region A the wave function will 
consist of solutions which decay as y — > — 00. These contain the terms 

exp(±^y a dy' p r y (y>fje X p(-~J V dy' p\{y')^ , (region A) (4.84) 
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y- y+ 

< B > 

FIG. 6 The classically forbidden regions, A, C, and the classically allowed 
region B for an excitation confined by A(y). 



where p l y (y) > and the phase reference point has been taken to be y_ 
Similarly in region C solutions will contain the terms 



exp 



i 



y+ 



dy' pUy') exp -- / dy' pUy') , (region C) (4.85) 



where now the phase reference point is Finally, in region B, the eigen- 
function will consist of a superposition of all four solutions corresponding 
to the four actions Sq (y) ( 4.82 ) , one for each of the classical momentum 
branches displayed previously in figure [|. 

Let us write the appropriate superposition in each region as "f^ (phase 
reference y_), iffg, ^ b (phase reference y+), and ^c- Our asymptotic 
solutions are only valid away from the turning points at y±. We match 
v[r c — > \p B and — ► \lr^ by analytically continuing our solutions into the 
complex coordinate plane in order to circumvent y± . The interested reader 
is refered to Appendix ^ for the technical details. Here we summarise 
the results. Starting from a complex evanescent solution satisfying the 
boundary conditions at y = +oo 



D 



(dE_\ 

\d P y ) 



u- QtI {y)e+W \ -±flpr y (y>)dy'-± S V y+ PW) d V 

voAy)e-^ 2 



Py iv) 

(4.86) 

we follow a contour in the upper half of the complex plane around y + to 
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obtain 



u+ 7 (y)e+^/ 2 

</(y)e-^ /2 



Py (a) 



Vz(y)e +#/2 

V(!/)e-^ 



Pv iv) 



(4.87) 



for y_ <y <y +1 and continuing around y_ we find 



%/(2/)e" #/2 A 1 



dp y 

i rv- 



p v (y) 



L4 



(sPy ) 



%x(») e_< * /2 



x e 



Py (y) 



£0 



x e 



OB 
9Py 

ft 



u+ I (y)e+^ 



Py (y) 

jy-pWW 



A 



dE 

<>l>!i 



<Av)^ /2 \(-i! y v + _Ptiy')dy' e -H:iPyWv' 
vUyy-^r 



Pv (y) 



x e 



-I r-pl(y'W +k S v y -pW)dy' 



(4.88) 



for y < y (The spinor and square root prefactors in ^ A and which 

are complex can be found in appendix |c], equation ( C.188| ).) To satisfy the 
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boundary condition at y = — oo that our solution vanishes, the coefficients 
of the first and fourth terms of ^^(y) must be zero, i.e. 

es/.-pjW dy'-ily- Py(v') d v' = e ^(2«+i)_ 



Our quantisation condition is then 



\ f + _ Py~(y') - Py{v') dy' = 2tt (n + , (4.89) 



which is precisely the quantisation rule used earlier (equation ( 4.72 )) to 
derive the spectrum of Andreev quasiparticles. We have also proved that 
the Maslov index is m = 2 (valid for quasiparticles bound by a smoothly 
varying order parameter, |A(y)|). 



5. THE SINGLE VORTEX 



The single vortex is an ideal problem to demonstrate how our semiclas- 
sical theory works since, as will be seen below, both the topological phase 
and the ft-dependence of the Hamiltonian play important roles. We will 
also (in an accompanying appendix) be able to make a direct comparison 
between the two semiclassical theories. 

The order parameter for a single vortex in a cylindrical coordinate sys- 
tem takes the form 

A(r) = |A(r)|e- ie , 

where the radial profile |A(r)| is shown in FIG. 0. Since [3] 

eA(r) B 

— ~ 1 

m/4> b c2 v ' 

for r < £, in this regime 



mv s (r) = 




i.e. the role of the superfluid is entirely represented through this ^-dependent 
term. Notice, since A(r) = the classical dynamics governed by Eg (po,r) 
is not influenced by the superfluid flow, whose role appears in a non-zero 
phase S"' J (r). However the classical dynamics governed by E a (p,r) does 
include the influence of v s (r). It is the latter which we investigate here. 

As already dicussed, HV(f)/2e can be interpreted as the vector potential 
associated with a line of phase singularities which in the present case runs 
along the z-axis carrying a flux, $o, from positive to negative z. From 
this vantage point we now investigate the behaviour of quasiparticles in 
the presence of a vortex. 
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|A(r) | 




FIG. 7 Profile of radial dependence of A(r) 



We follow Caroli, de Gennes and Matricon |29| by seeking a solution 
to the positive angular momentum branch of the spectrum, i.e. we take 
pe > 0. Since then only E + (p,r) > our Hamiltonian describing the 
quasiparticle-hole excitations is 



E + (Pr,Pe,Pz,r) 



Since it does not depend upon or z, pe and p z are cyclic variables (con- 
stants of the motion). p z is a continuous variable but pg must be quantised. 




5.1. pe quantisation 



We apply the generalised EBK quantisation rule ( 3.59Q to a path en 



circling the origin. For this path the Maslov index is zero, due to the fact 
that there are no turning points, but the topological phase associated with 
the vortex ( |3. 58 ) is important: 



1 ^ .„ _ / 



hj ^ = 2 M* 2 



where v and are integers. Writing po — hfi, /i an angular momentum 
quantum number, and using = w- § V</> • dr = —1, we find 



H = u+^, u = 0,1,2,. 
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The topological phase forces /i to be half integer. This agrees with the 
classic work of Caroli, de Gennes and Matricon pj| who obtained the 
result in a very different way, namely by matching wave functions. 

From our present vantage point we see that the effect of the half-flux 
flowing along the z-axis string is to give both particles and holes an angular 
momentum 'kick'. 

5.2. Radial behaviour of the excitation orbits 

Setting E = E + {p r ,pg,p z ,r) and inverting, the radial momemtum 
branches are found to be 



\ 



h 2 (fi 2 + i/4) 



± 2m 




the ± corresponding to particle-like or hole-like excitations respectively. 
Notice the effect upon pf{r) of having in the Hamiltonian is to 

introduced the two new terms fi- 2 /4r 2 and h 2 /i/r 2 . 

FIG. H shows the radial form of the orbit. One should notice the fol- 




FIG. 8 p^r{r) as a function of r. Solid line denotes a particle-like branch, 
a broken line a hole-like branch. The inset shows an enlarged view of the 
turning point where Andreev reflection occurs. 



lowing features: (i) The complete orbit is comprised of both particle- like 
and hole-like segments, (ii) As r — * both pf (r) — > i.e. we have classical 
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turning points. Since, by Hamiltons equations, 



= y '(E + h 2 n/2mr 2 ) 2 - |A(r)| 2 p ±( r ) 
E to ' 

for pf(r) > we have f + > and r~ < i.e. upon the positive quasiparti- 
cle branch the excitation is moving away from the vortex core whilst along 
the positive quasihole branch it is moving towards the vortex core, (iii) 
As the excitation propagates into the superconductor the velocity becomes 
zero at a point, r<i, defined by 

y/(E + h^/2mr 2 d f - \A(r d )\ 2 = 0, (5.92) 

and beyond this point p^r become complex. Thus is also a turning point. 
However since p^{rd) — p~(rd), rd is the point at which a particle-like 
excitation converts smoothly into a hole- like excitation (see inset in Fig. ||). 

5.3. Calculating the semiclassical spectrum 

Our general quantisation rule for the radial momentum reads: 

<j)pP(r)dr = 2nfr(n + , (5.93) 

where = ±. p@(r) stands for the appropriate branches of the momentum 
along the orbit in Fig. ||. (The Maslov index is to = 2 since the orbit is 
topologically a circle.) In general we cannot obtain an analytic solution to 



the integral in equation (5.93). None the less we can make some headway 
by using a simple model for the profile of A which retains the essential 
physics. 

5.3. 1. Model step profile for |A(r)| 

The model we adopt replaces the profile shown in Fig. |t] with a step 
profile. In particular we take: 

|AWI = {a! lit < 5 - 94 » 

where £ is a length scale characterising the distance over which |A(r)| rises, 
i.e. represents the vortex core. For instance we might take it to be the BCS 
coherence length £o = hvp /tvA^, but we need not do this. (We have in 
mind here the work of Gygi and Schliiter ||(| who have shown that £ -C £o 
as the temperature approaches zero, and indeed can end up comparable 
with the atomic spacing.) 
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Our model has the following features for the bound states E < : (i) 
when |A(r)| = Aoo the solution of equation ( 5.9 2| ) for the turning point is 



y /2m(A 00 - E) 



We note that r<i —> oo for /i — ► oo or E — > Aoo, i.e a bound excitation with 
in a high angular momentum state can in principle propagate far inside 
the superconductor, (ii) when |A(r)| = the only natural turning point 
encountered away from the core is set by £. Thus all excitations for which 
rd < £ are confined to the vortex core turning around at £, whilst those 
which have > £ penetrate into the bulk of the superconductor, turning 
around at r<i(fj,). 



5.3.2. Calculation of EBK integrals and Spectrum 

(i) For excitations confined to the core the p (r) simplify to 



Pt(r) 



, fi2(„ Tl /2)2 

Pf - Pi 2 ± 



(5.95) 



Integrals of these functions ca n be done exactly in terms of elementary 
functions. Thus the integral in ( 5.93 ) becomes 



\J {£,/L@f — 1 - arctan ^{^/Lp f-l 



where the length, L@, is defined to be 



p 2 F sin a + {32mE 



Here we introduced the angle a through p z — pp cos a. (a = or ir corre- 
sponds to the excitation travelling along the direction of the vortex axis.) 
Excluding excitations with high angular momentum or large p z (sin a — > 0) , 
we see this length scale is small compared to £. Thus, for (£/L^) 2 ^> 1, we 
expand our result to find 



pP(r)dr = Y,P 2h 



J /3 



7T 

2 



2e 



o((wc) 3 ) 



(5.96) 
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The energy, entering into this expression through L@, is quantised using 
equation (5.93). Taking the first two terms in (5.96), and seeking the 
energy in the form E = Eq + E\ -| — • we find for Eq 



ttHvfI sin a\ ( ml 
E ° = 2i [ n+ 4 2 



(5.97) 



Here we have used E <£ip 2 F sin 2 a /2m, so that our result is valid for states 
whose kinetic energy along the field direction is small. (This is consistent 
with the above approximations.) 

We are interested in bound states for which E < Aoo . To see how large 
Eq is we substitute into (5.97) with £o — tw>F /""Aoo and find (remember 
m = 2) 

_2 



E 



sin a n. 



(5.98) 



For states with sin a ~ 1 only the n = quantum number corresponds to 
a bound excitation. (This agrees with the findings of Bardeen, Kiimmel, 
Jacobs and Tewordt ||.) Notice our result becomes stronger for £ < £o, 
although for states with a large p z component new branches to the spectrum 
(n > 0) may appear. Thus to describe the low lying excitations in the 
vortex we must take n — and seek the next order correction to the energy, 
E\. Including the terms from ( 5.96 ), which should be evaluated at 

E = E = 0, we find 



2to£ 2 



In particular, for £ 



£o, the excitation spectrum becomes 

(?: 



E 



2 a 2 
A* — —■ 



(5.99) 



(5.100) 



We have recovered the famous Caroli, de Gennes, Matricon states |29| of 
low lying excitations confined to a vortex core which however were found 
by matching wave functions. 

As well as reproducing the excitation spectrum one can also show (see 
appendix |5|) that the radial behaviour of the semiclassical wave functions 
reproduces the correct r-behaviour near the origin where the exact result 
was also given by Caroli, de Gennes and Matricon Q . 

(ii) Now we turn our attention to states which penetrate into the bulk of 
the superconductor. For £ < r < rd we must add to our integral a further 
contribution with A(r) = Aoo. In this region p (r) are 



p k (r) = 



\ 



P 2 F 



h 2 {p? + 1/4) 



± 2m 




(5.101) 
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The analytic calculation of integrals of nested square root functions is very 
complicated. For example the Riemann surface upon which p + (r) (as de- 
fined by (§J0l|) ) is single valued has the topology of a sphere with 5 handles! 



The problem of evaluating a quantisation rule written as a line integral can 
be transformed into evaluating contour integrals on these Riemann sur- 
faces. In general, and in the example shown, the answer is not expressible 
in terms of either elementary or Elliptic functions. The integrals are how- 
ever rather easy to evaluate numerically. This we have done and quantised 
the areas to obtain the spectrum, E{p). The result is shown in figures [| 
and [h]. 



E/ | A (cd) 
1 r 



O.c 



0.6 - 



10 12 14 



FIG. 9 E{[i) flattens off as [i increases. 



Observe that whilst for small [i the spectrum increases steeply that 
once the excitation starts to penetrate the superconducting bulk it 'feels' 
less confined and hence E flattens off. Only very high angular momentum 
states have their energy approaching A^. 

As stated above, the small \i behaviour of E{p) agrees with the work of 
Caroli, de Gennes and Matricon |^|, and now we also see that the global 
behaviour of E{jj) is in agreement with the general findings of Bardeen, 
Kiimmel, Jacobs and Tewordt J^] . However our method of solution differs 
substantially from the work of these authors. 

6. SUMMARY AND DISCUSSION 

In this paper we have developed a semiclassical theory for quasiparti- 
cles in superconductors. In doing so we have pushed semiclassical methods 
to the limit by constructing a (fictitous) classical mechanics describing the 
orbits of quasiparticles propagating in the presence of lines of phase singu- 
larities. Adopting this approach enabled us to bring the full machinery of 
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FIG. 10 E(fj.) tends towards Aoo for very large angular momentum. 



modern semiclassics to the problem. In particular we used torus quantisa- 
tion to construct a generalised EBK quantisation rule for determining the 
scmiclassical spectrum. This rule included both the Maslov index, famil- 
iar from modern semiclassics, and a topological integer arising due to the 
global curvature of the space in which the quasiparticles propagate. The 
later is a direct consequence of phase singularities of the pairing potential, 
i.e., vortices. The power of this approach, first considered by Azbel' in the 
superconducting context, and extended by us here, lies in the general na- 
ture of the Hamiltonian system we have constructed. Unlike the approach 
of other authors [jl], ||, || where one must return to the BdG equations for 
each new problem and consider solving the differential equation across var- 
ious length scales, here our starting point is with the Hamiltonians, ( 3.41| ). 
and Hamilton's equations of motion. Once the pairing potential, vector 
potential and crystal lattice potential have been chosen for the problem at 
hand one can immediately investigate the quasiparticle orbits. In the case 
where the classical dynamics is integrable one can then proceed, via the 
generalised EBK rule, to quantise the quasiparticle motion. To make this 
possible a number of technical challenges needed to be overcome. Most 
important amoungst these was the observation that 'standard' semiclassics 
(i.e., asymptotics in the small parameter K) does not, in the multicompo- 
nent case, lead to the construction of a generalised EBK rule. We resolved 
this problem by constructing a new semiclassical theory whose correspond- 
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ing Hamiltonian system contains ft-dependent terms. For such a system, as 
we have shown, a generalised EBK quantisation rule does exist. We demon- 
strated the power of our approach by solving two well known problems cho- 
sen to elucidate those aspects of the theory which are new (^-dependent 
Hamiltonians and topological phases due to vortices). Of course the prob- 
lems we have considered have integrable classical dynamics. However there 
are many situations one may wish to solve where the classical dynamics 
is non-integrable. Thus we should examine which aspects of our theory 
remain valid in this case. 

When our Hamiltonian system exhibits chaos our corresponding (sta- 
tionary) Hamilton- Jacobi equation has no solution, i.e., the eigenstates of 
the BdG equations cannot have the 'simple' multicomponent WKB form. 
Since our 'classical' mechanics was derived in the first place using this 
ansatz for the wave function one must question whether the Hamiltonian 
system has any meaning for non-integrable systems. The answer is yes and 
the reason is as follows. If we start with the time-dependent BdG equa- 
tions, (replace E\ by ihd/dt in ( |2.l| )), and seek a solution in terms of a 
time-dependent multicomponent WKB ansatz, then we obtain, in place of 
the stationary Hamilton- Jacobi equation, the time dependent Hamilton- 
Jacobi equation: 

£°(VS(r,t),r) + — = 0. (6.102) 
A solution to this equation always exists unlike for the stationary version. 



The Hamiltonian system corresponding to (3.102) is identical to the one 
for the stationary equation and thus we have discovered our (fictitous) 
classical mechanics is the correct one to describe quasiparticle dynamics 
in superconductors regardless of the type of the dynamics exhibited. (The 
same cannot be said of our eigenfunction (^56|) , and the generalised EBK 



rule neither of which have any meaning in chaotic systems for which trace 
formulae must be derived fl.) 

One of our principle motivations for developing the scmiclassical theory 
was to construct a microscopic theory of quasiparticles in superconductors 
in large magnetic fields. Under such conditions the groundstate of a Type 
II superconductor adopts a so called Abrikosov flux lattice state j27|. 
This is a regular lattice of vortices, each carrying one superconducting 
flux quantum, $ , and having supercurrents associated with it. We have 
developed a simplified model of this state |5l| [32l [}3|, whose dynamics 
is integrable, and have used it to give an explaination of the origin of de 
Haas- van Alphen oscillations in Type II superconductors. 

Another problem for which the above strategy may be usefully deployed 
is that of the mesoscopic metals proximity coupled a superconductor. The 
refined semiclassics we have presented here may serve as a bases for in- 
cluding 'quantum diffraction' corrections missing from the naive analysis 
criticised by Altland, Simons and Taras-Semchuk jgjj]. 
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Finally we point out that the theory we have presented is restricted 
to s-wave pairing. An ongoing 'hot-topic' in superconductivity research is 
exotic superconductivity, in particular d- and p- wave pairing. We are in 
the process of extending our present theory to encompuss exotic pairing 
so that we may study quasiparticle orbits in the presence of d- and p-wave 
pairing potentials in magnetic fields. 

APPENDIX A: APPROPRIATE LABELS FOR THE ZEROTH ORDER 

AMPLITUDES 



Our zeroth order spinor amplitudes contain f3 = ± which are determined 
uniquely by j as follows: For a given branch Pg(r) of the multivalued 
momentum fuction we calculate the left hand side of 

Po , 1 



gL + _ mv ; + V (r) -e F = (3^{E? - p ■ vo)' - |A(r)|'. (A.103) 

Denote this quantity by T(pj). Then if T > 0, (3 = +, whilst for T < 0, 
= — . In this way each Pq determines a unique choice of (3. To give an 
interpretation to the sign of (3 we will need the following definition: If 

2 



> 



we call the excitation quasiparticle- like or simply 'particle-like', whilst for 

2 / . \ 2 



(«#(')) -(«oi(r) 



< 



we ca ll the excitation quasihole-like or simply 'hole-like'. From equation 
( |2.26 ) we have 

(<iw) 2 - (<i w) 2 - ^ Et -$:£: A * r ■ 

Noting sgn(E^ — p 3 Q -v ) — awe construct Table [l] to interpret the sign of (3. 
Since, as we have shown, /? is determined by Pq we can similarly interpret 
a given p^ as corresponding to 'particle-like' or 'hole-like' excitation. 

APPENDIX B: DERIVATION OF THE TRANSPORT EQUATION 
AND OTHER FIRST ORDER QUANTITIES 



Instead of expanding equation (2.11) directly we begin by taking its 
expectation value: 

( u* x (r) vl(r) ) x 

E x -H(p + ^,r) -|A(r)|e**W \/ u A (r) 
-|A(r)|e^W E x + H*{p+ ^,r) J \ v x (r) 

= F^DF = 0, (B.104) 
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Hamiltonian System 


Amplitude index 


Excitation type 


a = + 


/? = + 


particle- like 




P = - 


hole- like 


a = — 


P = + 


hole- like 




P = - 


particle- like 



TABLE 1 
Identifying the excitation type 



where we have introduced D for the matrix differential operator and F to 
represent the spinor. Then expanding the result, equation (B.104), upto 
and including first order in h we find 



= FjAjFo + FiD Q F Q + F^D F 1 + F$DiF . 



(B.105) 



Here the subscripts denote quantities of zeroth or first order in h. Dq is 
the zeroth order Hamiltonian matrix 



Do 



£i-tf e (p ,r) -|A(r)|e^W 
-|A(r)|e-^W E l + H$(p ,r) 



(B.106) 



i*o is the by now familiar spinor (2.9), D% is the first order matrix differ- 
ential operator, which we have not yet written down explicitly, but what 
is F\l Fi is the spinor obtained by going beyond the first two terms in 
the expansion of S(r) and S(r), equations and (2/7). Thus the wave 
function written to include the next order terms is 



ui(r) 

m(r) 



Uj.( r ) e iftS2(r) 
fJ I ( r ) e -iftS 2 (r) 

ui(r) 
vi(r) 



e zhS 2 (r) e *S (r)/h 



iS (r)/h _ 



«i(r) (+ift£ 2 (r)+z?iS 2 (r)) 
ik(r) (-iKE 2 (r) + ihS 2 (r)) 



o iS (r)/h 



(B.107) 



In (B.107) we used expfi^E^) w 1 + ihT, 2 , for iKS 2 <C 1 and so on. Turning 
our attention back to ( B.105| ) we recognise the first term as the zeroth order 
equation 

O (h°) = F$D F , (B.108) 
and because DqFq = the second term is also zero. This leaves 



o = f£d f 1 +f£d 1 f . 
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Since Dq = D$, the first term can, however, be written as 



F^D Fi = ( DIF ) F 1 = (D F y F, = 0. 



it; 



t 



Consequently the expansion of the matrix elements (13.104) gives the first 
order equation 



0{h) 



= F^D 1 F . 



(B.109) 



In particular notice that F\, the O (ft) correction to the wave function, does 
not feature in the expansion of the BdG equations up to first order in ft 
thus justifying the omission of such terms from the spinor elements ui(j) 
and vi(r) in (p79|) . 

To find the explicit form of D\ we require the order ft terms of H (p + ^f 1 
Writing H out we have 



in ^ dS ° 



V + P+ +V{r)-e F , 



1 

2m 



^(po,r) + -Lf-V-P^ 
2m V i 



ft_ 
P+ -V 

i 



(ft 2 ) 



where P+ = ^ + eA(r). H* (p + ^r, r) is written in similar fashion by 
replacing P + with P = ^M 2 - — eA(r). Then the matrix Di takes the form 



2m V % 



P+ . |v) 

+ i V*v-p- 



fv) 



so that equation (B.109) becomes 



FA 



i »VP+ 
o 



-^Mv-p- 

2 m ?. 



F 



F 



J_pH 

2m 







and if we then pull out the from the first term we obtain 



F = 0, 



-V • {FjM^} - - {vfJ} • M^ + ~F f M • V^ 



where the vector matrix M is given by 
M = 



= 0. 



+ f- 

2m 



(B.110) 



(B.lll) 
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The first term is purely imaginary since the differential operator acts on 
Fq Fq = \Fq\ 2 (M being diagonal). The remaining two terms taken together 
are purely real. To see this we rewrite the middle term in equation ( B.110| ): 



{V4} • MF = [(MF )t • (VF )] f , 



FjMt-VFo 



but since M = the last two terms in equation (B.110) take the form 



FqM. • VF 



= -2i Im 

i 



FqM. ■ VF 



which is certainly real. Equation ( B.llOj ) can therefore be separated into 
two equations which given explicitly are 



0(h) 
O(h) 




2Im( 



(B.112) 
(B.113) 



with = po(r)±eA(r). The first, (B.112), contains no phases and is the 
amplitude transport equation as we will show shortly. The other, equation 
( p3.113 ), determines the phase S[(r). 

In order to rewrite equations (B.112) and (B.113) in a more tractable 
form we will need Hamilton's equations ( 2.19Q 



dp 



and also the follwing relation 



( dEg(p,r) 
\ dr 



dE«(p,r) 
<9A(r) 



p=p?( r ). r=r °(*) 



(B.114) 



(B.115) 



where j Q (r) is the current density at the point r. We require an expres- 
sion for Hamilton's equations and the current in terms of the zeroth order 
quantities already found. For this we use the zeroth order equations (2.13) 



tf e ( Po ,r) |A(r)|e 



»0(r) 



|A(r) 



(r) 
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which we differentiate to give 



V dp 

( u?(r) 5f(r) )* 





uf (r) 
(r) 



(B.116) 



where 



= p -2Sl(r) 



(Ki) 2 + «l) 



2\ = p -2Sl(r) 



has been used. We can do the same for ^ncrv^ '■ 

9A(r) 



dA(r) 



= u 



0,1 



"fr ) 





(B.117) 



where the phases and amplitude, e 2Sl , have been cancelled because they 
shall not be needed. The explicit derivatives of H^ h are 



and 



dm 

dp 

dp 



dA 
dA 



1 P+ 

=— (p + eA = , 

m m 

= - — (p-eA) = - — , 

m m 



1 P+ 

= — (p + eA) e = e, 

m m 



1 



(p-eA) (-e) 



i-e), 



(B.118) 
(B.119) 

(B.120) 
(B.121) 



all evaluated for p = Pq (r). 

Now let us use Hamilton's equations and the current relation to rewrite 
t he tra nsport equat ion ( B.112| ) and the equation for ST (B.113 ). Using 
( |B .118 ) and ( B.119 ) we see that the right hand side of ( B.ll(j ) is p recisel y 
the inner product between the vectors in the order h equation ( B.112| ). 
Thus the transport equation can be rewritten as 



-2Sf{(r) UJ ^0 



dE%(p,r) 



dp 



P=Po ( r )/ 



Manipulating ( B.122 ) into the form 

' dE d(e- s i) 



E 



dpk dx k 



1 _ S i 



d 2 E 
dpkdx k 



0. 



= 0, 



(B.122) 
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we recognise the time independent transport equation of van Vleck |l4 
The solution is given by the determinant 



det 



drdl 



1/2 



(B.123) 



where c is a constant. 



All that is left to determine now is S{ (r). Returning to (B.113) we have 



= 2 Im ( v% 



ua ~ v a y( -^-(V5r + VEI)fif 
1 V +^-(TOr-VED«f 



o = {^ + (K I ) 2 -W 2 )£}-v^[ 

£ + (Ki) 2 - «i) 2 ) 



eA 

m 



V£[. (B.124) 



Now use equations (B.116) and (B.117), which written out in full become 



dEl 
dp 
1 dEg 
e OA 



eA 

m 



+ («i) 2 - Ki) 2 ) I 



m 

= - + (Ki) 2 - Ki) 2 ) — , 

m m 



(B.125) 



so that (B.124) becomes 

dEg 



dp 



V5[ 



OA 



or 



which we integrate along a trajectory of the Hamiltonian system from some 
initial time to to t 

Finally we have found an expression for S[: 

S[(v) = -- f j a (r).VSl< 
e ./to 



4G 



or 



1 f* , V4>(r) 



e 



S[(r)=-- r(r)--^dt'. (B.126) 



'to * 

We interpret this shift to first order in h of the phase in the following way. 



From (B.125) we have 

" \ r (r) = £ + (Kl)2 " W ' l)2) ^' (B ' 127) 

which is essentially the velocity of an excitation with an effective charge 
e*(r) = ((uqi) 2 — ( w oi) 2 ) e - This effective charge is determined by the 
relative amounts of particle and hole that the excitation is composed of. 
As the spinor is transported along the phase space trajectory, figure |l], the 
particle-hole composition slowly changes giving rise to a changing charge 
e*(r) (charge is not a good quantum number in superconductors) and hence 
a changing current j Q (r). This varying current interacts with the phase 
gradient of the order parameter and modifies the phase of the wave function 
to first order in h. Notice in particular that if no phase gradient exists, i.e. 
<f> = constant, then S\ = 0. 



APPENDIX C: TECHNICAL DETAILS OF THE ANALYTIC 
CONTINUATION OF THE SNS JUNCTION EIGENFUNCTIONS TO 
DETERMINE A QUANTISATION CONDITION 

The basic idea of the analytic continuation of solutions is simple: our 
asymptotic solutions are valid not only along the real axis but for complex 
coordinates as well, so long as we stay sufficiently far away from y±. Thus 
we obtain an approximate solution to the BdG equations throughout the 
complex plane (excluding discs, radius p say, centred upon y±) by contin- 
uing our asymptotic solutions to complex coordinates. In particular we 
can continue a solution valid for y > y + + p along a semicircular path, 
radius > p, in the plane and arrive at y < y + — p, thus determining the 
appropriate amplitudes and phases of the solution for y_ + p < y < y + — p. 
Similarly we can circumvent y_ to obtain the correctly phased solution for 
y < y_ — p. The situation is however complicated by Stokes phenomenon. 

Stokes phenomenon is named after its discoverer Sir George Gabriel 
Stokes (1819-1903). A general discussion of historical issues can be found 
in Heading j35|. Stokes, studying Airy's equation, noticed that if a general 
solution, given in terms of the two power scries with arbitrary coefficients, 
was represented by a certain linear combination of the asymptotic solutions 
in one region of the complex plane then in an adjacent region the same 
general solution was represented by a completely different form of linear 
combination of these asymptotic solutions. He discovered the constants of 
the linear combination changed discontinuously when crossing certain lines 
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in the plane - now known as Stokes lines. (Note it is a change in form not 
the numerical value which takes place.) 

For us, this means that if we choose an asymptotic solution, i.e. we set 
the coefficients of other solutions to zero, and then analytically continue, 
we will cross Stokes lines where the coefficients can jump, in particular 
from zero to non-zero, so that solutions which were absent suddenly appear 
changing the form of our solution. What we need to know is how to locate 
the Stokes lines and what rule determines the change in the coefficients as 
the lines are crossed. 

C.l. Stokes phenomenon in the presence of two solutions 

Stokes phenomenon occurs when the exponentially subdominant solu- 
tion out of a pair is at its smallest. Consider the asymptotic (WKB) form 
for the solutions to Airys equation: 

For real z these solutions are oscillatory but for complex z the solutions 
aquire the exponential factors e + ' tu ' z ^ or e^'" 1 ^-", with w{z) = |Rc{iz 3 / 2 }. 
Thus one solution is exponentially dominant, the other exponentially sub- 
dominant. Now e + \ w ( z ^ is maximally dominant when Im{iz 3 / 2 } = 0, and 
it is then, when the second solution is least visable, that its coefficient can 
jump. Thus Stokes lines are defined by 

lm{iz 3/2 } = 0, (C.128) 

in the present case. 

Anti-Stokes lines are defined by 

Re{iz 3 / 2 } = 0. (C.129) 

These are the lines along which the effect of the Stokes jump becomes 
important because both solutions are purely oscillatory. The solutions are 
said to be neutral, neither being subdominant wrt. the other. Notice that 
both the Stokes lines and anti-Stokes lines eminate from a branch point at 
z = 0. 

The general prescription when crossing a Stokes-line is given by 

[New subdominant coeff.] =[01d subdominant coeff.] 

+ (Stokes constant) x 
[coeff. of dominant term]. 

By following a given solution along a circuit surrounding the branch point 
a number of Stokes constants will be introduced. Requiring our intial and 
final solutions to be the same determines the value of these constants. 

The situation we must solve involves four solutions. Before we can 
proceed we must therefore generalise the above discussion. 
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C.2. Stokes phenomenon in the presence of more than two 

solutions 

[] Suppose now we have N solutions. Let us change our notation and 
write each solution as 



(C.130) 



i = 1,...,N. Then one might expect that as we move around in the 
complex plane a given solution, say i, could pick up contributions from 
some or all of the other solutions, j =/= i. The question we must ask is when 
can Stokes phenomenon occur? The answer is as follows. Introduce the 
singulant for a pair of solutions as the quantity <f> t — <pj . Then a necessary 
condition for ^ to experience Stokes phenomenon is 



lm{(f>i - <f>j} = 0, i^j, 
and sufficiency ^ is ensured by 

Re{&} > Re{0j}. 



(C.131) 



(C.132) 



We can understand the first of these conditions, ( |C.13l| ), as being a gen- 
eralisation of the Stokes lin es for t wo solutions. Indeed when we have two 
solutions so that <pj — —<f>i ( C.13l[) reduces to 



Im{2&} = 0, 



which is just the usual equation ( p. 128 ) for Stokes lines. When there are 
more than two solutions, so that in general cfij ^ 0;, (C.131) says that 
Stokes lines exist where the imaginary part of the exponents are the same. 
This ensures that when comparing two solutions the effect of the phase 
can be discounted. Then we can concentrate on the real parts which must 



satisfy (C.132). Clearly (C.132) implies 



so that is maximal over ifrj whose coefficient can then change. 

Thus Stokes phenomenon occurs where solutions are pair-wise maxi- 
mally dominant and subdominant. 

Anti-Stokes lines are defined by 



Re{& 



,}=0. 



(C.133) 



This situation will be clarified in the specific treatment of the semiclassical 
SNS junction eigenfunctions which follows. 

7 I am indebted to Chris Howls for explaining how Stokes phenomenon is to be un- 
derstood and treated in the case of more than two solutions. 

8 Assuming the Riemann sheets on which the individual solutions live are not locally 
disjoint. 
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C.3. Conventions and notation 

It will be useful to introduce the following compact notation for the 
semiclassical solutions with phase reference y + : 



\d Py ) 



Mz)= ^{Zti^) e v+ ' (c - 135) 



(^Ph ) 



P» («) 



Pb (z) 



where the &i(z) are analytic continuations of the ^i{y) to complex coor- 
dinates in a domain suitably cut to ensure single-valued definitions (see 
FIG. fj"l| and the next section). The same four solutions written with phase 
reference y- are obtained by replacing y+ by y_ in equations (C.134)- 
( |C. 137 ) and shall by distinguished from ^i{z) by appending a minus sign 
as a superscript thus 

*rw=- If (das) 

etc. Note that 

where we have introduced [1] as the pure phase factor 



We also have 



[i]=exp (U* 

[2] = [1]*, 

[3] = cxp (U 

[4] = [3]*, 



dz' p+(z') 



v+ \ 
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where * indicates complex conjugation. Furthermore we use <j>i and <p i 
(i = 1, . . . , 4) to indicate the exponents including the asymptotic parameter 
thus 

^= l -£ + dz'p+(z'), (C.139) 
= \f~ dz 1 p+(z'), (C.140) 

and so on. 



C.4. Ensuring Py(z) and the wave functions are single- valued 



In order to construct a solution valid in the complex plane we must 
have single- valued WKB solutions. The momentum branches for the su- 



perconducting state are multivalued functions. Recall pt (y) (4.67): 



P P v (y) = \jp F ~Pl-Pl+ P2m^E 2 - | Agp- 
In fact Py(y) are two branches of one function p y (z) given by 

Py(z) = \J P 2 F - P% - P 2 Z + 2my / E 2 - A 2 (z), 



(C.141) 



which is single-valued on its Riemann surface comprising more than one 
Riemann sheet. (Here we use A(z) to represent the analytic continuation of 
|A(y)| to complex values.) Suppose we have the branch of (C.141) which 
is equal to Py(y) on the real axis, then we obtain Py(y) by analytically 
continuing the function around one of the zeros of \J E 2 — |A(y)| 2 . We can 
see this explicitly if we consider a smooth slowly varying |A(y)| and expand 
in the vacinity of the turning point 



|A(y)| = A 



9|A| 



dy 



(y-y±) + 



(C.142) 



II ± 



The function ^JE 2 — A 2 (z) then becomes (A = E) 



VE 2 -A 2 (z) = J2E\W y± Q± e 



(C.143) 



where z± = g±e ie± t 2 are the coordinates centred on y±. Clearly 6± 



9±+2tt causes the function in ( p. 143 ) to change sign so that p y (z±) = Py(y) 



becomes p y {z±e l2lT ) = Py(y)- To make p y (z) single- valued in the complex 
plane we insert a branch cut between y+ and y_. Notice that, so long as 
Pp—Px—pl > 2mRe{y/E 2 — A 2 (z)}, the Riemann surface of p y (z) consists 
of just two sheets. However our WKB solutions are constructed from p^ 
and — Py so that we must consider the four Riemann sheets of p y (z) and 
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z-plane 




FIG. 11 Branch cuts which ensure the WKB solutions are single- valued. 



—Py(z) together. FIG. [TT] shows the plane in which the WKB solutions 
are single- valued. The second branch cut emanating from y_ is necessary 



because from (4.81) 



f dE 

\d Py 



which in the vicinity of y_ is 



-1/2 



ys 2 - a 2 (z)' 



(C.144) 



BE 
dp y 



-1/2 



-iB-/i 



*y/2E\*£\ V _ Q. 



(C.145) 



Now (C.14E) multiplies each of the WKB solutions and changes sign when 
6 — > 9 + 4:7: even though the momentum is single- valued under this change. 
The WKB solutions are then only single-valued if a second branch cut is 
inserted (as shown). 



C.5. Rules for continuing across branch cuts 

Let us start with the turning point y+ (FIG. |Ti"| ), and consider p y (z) — 
Py{z). By inserting the branch cut between y+ and j/_ we have prevented 
p y {z) taking the value Py(z) at each z. To display all of the values of 
p y {z) we take two copies of the complex plane, label them sheet 1 and 
sheet 2 respectively, and assign sheet 1 as the domain of p+ and sheet 2 
the domain for p~. Look at Fig. [l^. If a contour respects the branch cut, 
71 for instance, it stays on the first sheet. However if we follow a contour 
like 72 when we cross the branch cut of sheet 1 p y {z) is changing smoothly 
and so we appear on sheet 2 and p y (z) has become Py{z). What happens 
to the eigenf unctions? 
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0=271 






0=0 







FIG. 12 Alternative contours in complex plane. 72 crosses the branch cut. 



Consider 
*i(z) = 



where 



f dE 

\dpy 



-1/2 



#1/2 



and the spinor is given by (see equations (4.79)-( p~80| )) 

/ 



V 



(C.146) 



(C.147) 



(C.148) 



f3 = + corresponds to p+ . Now suppose we follow the two contours shown 
in FIG. The first, 71, does not cross the branch cut. Let us use z a bovc = 
re i0 a bovc t label the z-coordinatc in the plane from following this contour. 
If instead we follow 72 we do cross the branch cut. Let the coordinate from 
following this path be Zbeiow We then have z a bove = Zbciowe +l27r . What 
we require is 'l'i(zbeiow) written in terms of z a t>ove- Let us consider in turn 
each of the terms entering into ^i(z). 

C.5.1. The change in (dE/dp y ) 

Firstly consider ( |C.147 ) and note that 



V^ 2 - A 2 (z) 



2-E| % y \y + 



r e 



i0/2 



(C.149) 
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in the vacinity of y+ so that ( |C.147 ) becomes (using Zbdow) 



dE\- 1/2 ^l/2 e -^bo!ow/4 ! 

— ■ ■ (C.150 



We then have 



fdE_y 1/2 fdE_y 1/2 



\ d Py J p + ( 2bolow ) V d Py J P + (z abovo e-* 2 " ) 

£l/2 e -i(0 aboV e-27r)/4 



W 2E \%t\y+ r ^Jp+ (z ahovc e-^) 

= ein/2 (^-) 1/2 ' ( C - 151 ) 

where Py(z a bovce~ l2lr ) = Py (z a bovc) has been used. We find this amplitude 
aquires a factor +i. 

C.5.2. Change in the spinor 

Here we also have the spinor ( C.148Q which using ( p. 149 ) becomes 



< / (zbciow)e+^ /2 \_f ^/(^abovc)e+^/ 2 . 
< J (^below)e- i ^ 2 J " ^ %/(^bove)e-^ 2 J ' 

C.5.3. Change in the phase integral 
Finally consider the phase 

i f v+ 

</>l(zbclow) =T / Py(z')dz', 



' ^belo 



p+(z')dz'. (C.153) 



Changing the limit z a bovce _l2,r to z a bovo will mean the phase of z' along a 
contour is decreased by 2ir, we must therefore explicitly include e~ l2jr to 
maintain the value of ( p. 153 ) thus 



<M^bciow) = \ I p+(z"e- i2 *)dz", 



n 



Zabo 



p-{z")dz' 



). (C.154) 
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C.5.4- How the eigenfunctions change 



Putting (C.151), (C.152) and (C.154) together we find the rule for cross- 



ing the branch cut to be 



( Mi\ ^ v oj( Zaho 

\ d Py) _, , 

Py (Zabovc) 

-«* 3 (2abovo)- (C.155) 



Clearly we also have 



*2(2bolow) = +«*4(^abovc), (C.156) 
I^Obclow) = +«*l(Zabovo), (C.157) 
* 4 (2bclow) = +i* 2 (2abovo)- (C.158) 

These rules apply when crossing the branch cut in the clockwise sense 
around y + . When moving in the positive sense +i is replaced by —i. It is 
also easy to see that ^j" = etc. for solutions with phase reference 

y— 

C.5.5. A rule for the pure phase factors 

We have previously introduced the definition 

[1] = cxp U J"' dz' p+iz'^J . 

The integral may be taken above or below the cut connecting y+ to y_. 
Since p+(z abo vc) = Py(z bc i ow ) we have 

[ljabovc — [3]bclow> 

or 

[1] - [3]. (C.159) 

Similarly: 

[2] - [4], (C.160) 

[3] -> [1], (C.161) 

[4] -> [2]. (C.162) 

We now have all the rules needed for crossing the branch cut. 

Our task is now to find the location of the Stokes and anti-Stokes lines 
as defined by fl013l| ) and ( |CT33| ). 
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6=271 



9=71 



e=o y + 

FIG. 13 Branch cut in complex plane. 



-*■ y 



C.6. Locating the Stokes and anti-Stokes lines 

In order to locate the Stokes and anti-Stokes lines in the vicinity of y+ 
and y_ we require the form of the momentum branches there. These are 



\ 



P' r -Pl-Pi±2m. IE 





dA 




2E 


z 




dy 


y± 



(C.163) 



where we have used the expansion for A(z) (equation ( C.142| )), and z ~ re l6 
is the coordinate in the complex plane centred on either y + or y_ (which of 
these it is should be clear from the context). Firstly consider the vicinity 
of y+, FIG. [l| We require the exponents of the WKB solutions. Let us 
start by considering 4>\ 



y+ 



p+(z')dz', 



v+- 



p+(z)dz. 



(C.164) 



To calcul ate the form of this integral we note that p%— p 2 — p 2 3> 2my / 2E\ If |i 
so tha t (|C.163| ) can be further approximated by expanding in powers of 
z 1 / 2 . ( |C.163 ) then becomes 



Py{z) = /' 



dpi 



dz 1 / 2 



' 2J5| |f I 



(C.165) 



Only the 



where p^ — y 'p 2 F — p 2 — p 2 and p^ = 

first two terms of ( C.165| ) shall be needed. Inserting the p+ expansion into 
( |CT64| ) we have 



MZ)= h 



y+-z 



or using z 



dz^+p^z 1 ^ , 

= \ {p {0) (»+ - *) + \p {1) (y+ - z) 3 - 

= along bottom of the branch cut in FIG. |l3|) 



n i o 



(C.166) 
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If we write 

{}± = {p< >re M ± H^V^e**/ 2 } , (C.167) 

then the four exponents of the WKB solutions with phase reference y + can 
be written compactly as 





(C.168) 


2(Z) = -{{} + , 


(C.169) 


.W = +j{}-, 


(C.170) 


4^) = ~0- 


(C.171) 



C.6.1. Stokes and anti-Stokes lines around y + 

For the Stokes lines we require Im{</>j — cf>j} = for each distinct pairing 
of ( C.168 )-( C171 ). We consider one example and then the remaining pairs 
have their results summarised in table ||. 

Consider <f>\ — <j>2'- 

i — i 

- w 

The imaginary part is then 



Im{<f> 1 -M = l(p i ° ) rcos8 + ^pWr 3 / 2 cos(?f) ) . (C.172) 
whilst the real part is 



2 



Re{& -&} = -| (p {0) rBme+lp^ 2 sm(jf) J . (C.17>!) 



The Stokes lines are determined by the angles 9 for which Im{</>i — 02} 
vanishes. It will be sufficient to locate these lines approximately. Close to 
y + (i.e. r — > 0) we can in the first approximation neglect the second term 
in flC.1721 ). Then 

p (0 Vcos6» = 0, 



gives 0° = vr/2,37r/2. Returning to flC.172| ) and substituting 9 = 9° + 69 
we can determine how the second term shifts 9. Since 69 is small all we 
require is its sign. For 9° — n/2 we find 69 = — \S9\ and the Stokes line is 
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FIG. 14 Stokes and anti-Stokes lines around y+. 



then located at 9 = n/2 - \56\. For 9° = 3ir/2 we also have 59 = -\60\ as 
can be checked. Thus 9 = 3tt/2 - \S0\. 

To decide which function \fi or ^2 is dominant and which is subdom- 
inant on the Stokes line we consider the sign of Re{0i — 02 }• For this 
purpose it suffices to ignore the second term in QC.173| ). For 9 = tt/2 - \66\, 
sin(vr/2 - \S9\) > 0, Re{0i - 2 } < so that * 2 > At 9 = 3tt/2 - \S9\ 
the opposite is true > ^2)- 

We have now located the Stokes lines for the pair of functions ^1 and 
^2, and we know the dominancy on each line. Calculating the anti-Stokes 
lines proceeds similarly. The results are summerised in table ||. 

FIG. [li] shows the Stokes and anti-Stokes lines around y+. Notice that 
the anti-Stokes lines at +tt — \S9\ and + tt + \S 9\ are shifted off the real 
axis (\59\ =^ 0) by the second term in ( |C. 173 ) . These lines are where 
4>i — 2 and 4>3 — 04 are neutral but since 2 — — 4>i, and 04 = ~03 
this neutrality condition reduces to ^ 1,^2,^ 3 and ^4 being individually 
neutral (i.e., oscillatory) on the appropriate Stokes line. It is these Stokes 
lines no longer coinciding with the real y-axis that ensures the existence of 
evanescent solutions for y — > +00, i.e on the real axis. 
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TABLE 2: Summary of Stokes lines and anti-Stokes lines located around the turning point y + . 



C.6.2. Stokes and anti-Stokes lines around y_ 

So far we have only considered the turning point at y + . The situation 
is very similar at j/_ . Note however that the phases are given by 



V 

x—V- 







Py(z)dz, 



-Up® { Z -y_) + 2 ~pV (z-y.f 2 



(C.174) 



where now z — y- = re l ° (9 = below the cut again), i.e. 



;{}+. 



(C.175) 



and 



-(*) = + 



{}-, 
{}-, 



(C.176) 
(C.177) 
(C.178) 



(compare (C.16S)-(C.171)). Using these the Stokes and anti-Stokes lines 
around y_ can be calculated. The results are summarised in table and 
drawn in FIG. The z-plane containing both y + and ?/__, and all 

the Stokes-lines and dominancy changes is shown in Fig. [H| This figure 
contains no fewer than 18 Stokes lines and 11 anti-Stokes lines. When 
there are more than two solutions it becomes difficult to label a solution 
with dominance or subdominance. Dominant with respect to which of the 
other solutions? The approach we have taken here is to write on the figure 
either > vE^- (where we use the > symbol not just to represent greater 
than but rather to indicate maximal dominance of ^ over ^j) or (\E^, ^j) 
which indicates that the pair of solutions are neutral with respect to each 
other. Clearly > tyj is needed at each Stokes line and (^j, at each 
anti-Stokes line. Thus if a > <J/2 Stokes line is crossed the coefficient 
of $2 undergoes Stokes phenomenon unless is not present. Notice that 
along the real axis for y > y + and y < y^ we have coincident Stokes and 
anti-Stokes lines. At first this seems like a contradiction but it is not for 
the following reason. The anti-Stokes line means that the pair (^1,^4) 
is neutral, as is (\&2, ^3)- But according to the Stokes line ^1 > ^3 and 
* 4 > * 2 - Taken together these statements imply (^1,^4) > (^2, ^3), 
or in words, the neutral pair ^4) is maximally dominant and (^21 ^3) 
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FIG. 15 Stokes and anti-Stokes lines around y_ 
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FIG. 16 Location of the Stokes and anti-Stokes lines for the SNS junction 
problem when there are four WKB solutions to consider. 
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TABLE 3: Summary of Stokes lines and anti-Stokes lines located around the turning point y_ 



maximally subdominant. Explicitly the solution along y > y + which decays 
is 

* w>1/+ (y)=A*2 + B*3, (C179) 



(JUL) 



v+j(y)e-^ 2 



+ 



Py{v) 

B ( ii 



Vor(y)e- t4>/2 



Py (v) 



(C.180) 



where the modulus in the exponent ensures an evanescent solution. In 
( p. 180 ) we have used p r y and p y to represent the real and imaginary parts 



respectively of the momentum. Note that p y = (p y )* so that the real parts 
are the same. Now we can understand fully the more general definitions 
of the Stokes lines and anti-Stokes lines. The pair (^2, ^3) are neutral 
because they have a common damping factor, neither is more damped 
than the other. It is only the real part of the exponents, <f>i,4>j, which 
control this hence Ke{cf>i — <f>j} — is the general anti-Stokes condition. 
The exponentially growing solution is 



^y>y + (y) = C^l+D^ 4 



(C.181) 



c (u+Me +i4,2 \ -UIpWW +h Iy + pl(y'W 

dPv) Pi(y) 

d u^ +iH2 \ e HIlpWW e H Iy + pi(y'W 

Vor(y)e- l4,/2 



+ 



(MR) 

\d P y J 



Py (y) 



(C.182) 



Comparing this with ( C.180| ) we see that <3>i and 'J 3 have the same oscilla- 
tory part, likewise ^4 and ^2- Again attention is focused on the real part, 
Re{(/)i i 4} > Re{(/>2,3}, because we can discount the common phase. The 
condition for discounting the phase is the general Stokes-line statement 
Im{^ - fa} = 0. 

Suppose we start along y > y + with the decaying solution ( C.179| ). If 
we follow the solution around y + in the upper half plane (see FIG. |14| page 
p8| ) we expect the first Stokes jump to occur at 9 = +37r/2. Continuing 
around to 9 — +2tt we expect to have had four Stokes jumps and hence 
four Stokes constants are introduced. All four solutions are then present 
on the real y-axis between the turning points. The solution along y < ?/_ 
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FIG. 17 Stokes and anti-Stokes lines around y+. 



would in principle contain 10 Stokes constants. Fortunately there is a 
quite remarkable simplification as will be seen in the next section when we 
calculate the Stokes constants. 

This concludes our discussion about Stokes and anti-Stokes lines. 



Stokes constants can be calculated by following the changes of a solution 
when passing all the way around a turning point and back to the starting 
point. 

We will now calculate the Stokes constants corresponding to the Stokes 
lines in figures [l4| and [l5| We shall follow the changes in a solution, evanes- 
cent along y > y + , as we move around y + in the complex plane. We will 
need to know what to do when crossing the branch cut along y < y + . 

FIG. ^| shows the various sectors in the complex plane, and we have 
also included a Stokes constant corresponding to each condition "J,; > \&j 
which is relevant when starting from the evanescent pair (\&2, ^3)- Notice 
that while we require the solution for y > j/+ to decay in order to satisfy 
the physical boundary conditions any arbitrary combination of ^2 and ^3 
have this property. Thus we consider = A^2 + B 1 ^^ (A, B arbitrary). In 
what follows the numbers indicate the sector under consideration and next 



C.7. 



Calculation of Stokes constants 
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to them we have written the form of the solution in that sector. We have: 



1 - 4 : A3> 2 + -6*3 

5: (i+BQ)f 2 +M 3 

6: (A + BQ)^2+ BV 3 +Bm 4 

7 : SBV 1 + (A + BQ)^ 2 + B^ 3 

+ (BR+(A + BQ)T)^ 4 
8: -iSBV 3 -i(A + BQ)t> 4 -iB^ 1 

-i(BR+(A + BQ)T)^ 2 
9: -iB(l + SU)^x + same * 2 ,*3, 

-i(A + BQ + (BR + (A + BQ)T)V)^i 
10: -i(B(l + SU) + W(BR+ (A + BQ)T))^ X 

+samet 2 ,*3,*4 (C.183) 
11: -i(B(l + SU) + W(BR+(A + BQ)T) 
+Y(A + BQ + (BR +(A + BQ)T)F))*i 
+same ^2,^4, 

-i(SB + X(BR +(A + BQ)T))^ 3 
12-14: -i(B(l + SU) + W(BR+ (A + BQ)T) 

+Y(A + BQ+ (BR +(A + BQ)T)V))^x 
-i(BR+(A + BQ)T)V 2 
-i(SB + X(BR +(A + BQ)T) 

+Z(A + BQ+ (BR +(A + BQ)T)V))V 3 
-i(A + BQ + (BR +(A + BQ)T)V)V 4 

Comparing the solution in sector 14 with that in 1 fixes the Stokes con- 
stants. To help us here we are allowed to treat each equation obtained by 
equating coefficients as two equations because the parts depending upon A 
and B can be matched separately (remember A and B were arbitrary so we 
can vary either independently) . Thus for instance matching the coefficient 
of *2 gives 



or 



A = -i(BR+(A + BQ)T) 



1 = — iT (from varying A), 
= R + QT (from varying B). 



(C.184) 
(C.185) 



( |C . 1 84 ) gives T = +i. Proceeding similarly we find 



S = T = U = V = +i, 
Q = R = W = X = 0. 

Y and Z are undetermined but this does not matter because we can avoid 
crossing the corresponding Stokes lines when solving the two turning point 
problem. 
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FIG. 18 Stokes and anti-Stokes lines around y_. 



Calculating Stokes constants around y_ proceeds in exactly the same 
way. FIG. [l^ shows the various constants. (The extra branch cut between 
sectors 14 and 15 means *~(zi 4 ) = — ^~(zi 5 ).) Note however that moving 
through the sectors 1 — > 14 we are moving around j/_ in the negative 
(clockwise) sense and we require Stokes constants in the positive sense. 
Taking this into account we find 

E = F = G = H = +i, 
C = D = I = J = 0. 

K and L are undetermined, but again they will not be needed. 



C.8. Discussion of Stokes constants 



We have found that all the constants are zero apart from those corre- 
sponding to 9 — +7r/3,+57r/3 around y + and 6 = —ir/3,+7ir/3 around 
y_. This is surprising since it appeared that the sufficient condition for 
Stokes phenomenon to occur was given by ( |C.132 ). If we consider again 
going around y+ we find the following is true 

^3 only 'sees' ^i, 
^2 only 'sees' ^4. 



Return to the definitions of Vf^, equations (C.134)-(C137), then we observe 
that the exponent of Vti depends upon +Py , and \&3 upon +p~ , whilst '5 2 
and ^4 depend upon — p+,and — p~ respectively. This is the resolution of 
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FIG. 19 Location of Stokes lines and anti-Stokes lines once the simplifica- 
tion due to disjoint Riemann sheets has been taken into account. 



the problem. When analytically continuing around the complex plane Py~ 
can be continued to p~ and likewise — Py to — p~ , but there is no way (at 
least locally in the plane) for p+ to be continued to — Py . This would require 
the real part of Py to pass through zero, but we have assumed p 2 F —p\—p z 
is large so that ±2m^jE 2 — A 2 (z) cannot reduce the momentum to zero. 
However one might imagine that |A(y)| analytically continued into the 
complex plane may have singularities and might therefore reduce Py to 
zero somewhere. For the moment at least we had better say that 

if locally in the complex plane solutions cannot be analytically 
continued one into another then there will be no Stokes phe- 
nomenon between them. 

To put it another way, solutions on disjoint pieces of a Riemann surface 
do not experience Stokes phenomenon. Note that p y {z) has two sheets 
(corresponding to p y (z = y) = Py{y) and p y (z = y) = Py(y)) and —p y (z) 
also has two sheets — p+ and — p~ , but although all four sheets make up 
the Riemann surface for the solution the sheets exist in mutually disjoint 
pairs. 

C.9. Derivation of a generalised Bohr-Sommerfeld quantisation 

rule 

Returning to the problem at hand it is clear that a significant simpli- 
fication of FIG. |l6| is possible. We need only consider those Stokes lines 
(shown in figure |l9| ) whose Stokes constants are non-zero. These Stokes 
constants are all +i when the lines are crossed in the positive sense. Fi- 
nally then, let us follow the evolution of an evanescent solution constructed 
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out of ^2 and ^3 through sectors 1 — * 6. We have: 

1 - 2 : A^ 2 + 5*3 

3 : iB^x + A^ 2 + B^ 3 + iA^ 4 

4 : iB[l]«7 + A[2]^ + S[3]*3 + iA[4}^ 
5-6: iB ([1] + [3]) *r + A[2]*a + B[3}% 

+ »A([4] + [2])*4. (C.186) 

Now for y < y- the evanescent solution must consist of Vf^ and *g only 
so the coefficients of and $4 must both be zero. For this yields 

[1] + [3] = 0, 
[3] 



or 

g h Jy- I'll""' Uj-fjH'"» = g 

We have derived the quantisation condition 



%Iy-Pt(v') dy'-j-f^p-(y') dy' = i*(2n+l) 



l j V+ _Pi{y r )~Py{y') dy 1 = 2n(n+iy (C.187) 



used earlier (equation (4.72)) to derive the Andreev spectrum but we have 
also found the turning point correction 7 = 1/2, i.e. the Maslov index 
m = 2. This is what is expected for the Maslov index of a Lagrangian 
manifold which is topologically a circle. In the present case, from our 
approach, we see each turning point introd uces a phase change of 7r/2, 



i.e. contributes 1 to the Maslov index (see ( C.187 )), whenever the local 
topology of the Riemann sheets in the vicinity of the turning point consists 
of two sheets. This concludes our derivation of the quantisation condition. 

Before leaving this section we wish to comment that in the main body 
of this paper the explicit solutions are written down for each region. Since 
we require these for z = y, between the turning points, and our solutions in 
sector 3 and 4 are for z = ye t27r , we must take the functional forms ^i(?/e l27r ) 
and rewrite them in terms of ^i(y) i.e., we must follow the solutions in 
sector 3 across the branch cut connecting y + and Again we use the 



rules flC.155| )-( |CT58l), and ( |C.159| )-( |C~1~62] ). Thus for example the solution 



^b(u)) equation (4.87), was obtained by following the solution in sector 3: 

*fl(«above) = IB*! + + B^ 3 + iAV 4 , 

across the cut to obtain 



08 



and then explicitly evaluating this at Zbciow = V- 
The prefactors 

dE_y 1/2 ( 4Me«+'* \ 

in the regions y < y_ and y + < y are complex and can be shown to be 

(dEY 1/2 ( ulMe^ 2 ) _ a(y) ( e+^+^ \ wu-w 
UJ P g { <Ay)e-^ ) - a(y) { e-WM-W ) 6 



(C.188) 



where 



if mV2 \f E 2 V 1/4 f \A(y)\ 2 - E 2 \- 1/S 



<%) = ^arctan — — . 



APPENDIX D: LIMITING BEHAVIOUR OF BOTH SEMICLASSICAL 
WAVE FUNCTIONS AND THE EXACT SOLUTION AT THE ORIGIN 

OF A VORTEX 

We now show that the limiting, r — > 0, b ehaviour of both the semiclas- 



sical wave function presented in section 2JB an d th e wave function of the 
effective semiclassical theory derived in section 3.2 agree with the known 
exact solution of the Bogoliubov-de Gennes equation at the origin of a 
s-wave vortex. 

D.l. Asymptotic behaviour of the exact solution 

It is well known jl| that in the r — > limit the BdG equations can be 
solved exactly. Writing the wave function in the form 



the BdG equations are reduced to a radial equation for /(r): 
-h 2 (l d d If 1\ 2 2mE S 



°z—\~J^~-A»-°>o) + M=0, (D.189) 
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(since |A(r)| — > as r — > 0.) Here k 2 = k 2 F — k 2 , kp the Fermi wave vector, 
and <7 Z is a Pauli spin matrix. This equation has the J v (z) Bessel functions 
as solutions: 



f(r) 



A+Jft-i/2 ((hp + q)r) 

A-J^+i/2 {{hp - q)r) 



(D.190) 



where, /i ± 1/2 are integers, and following de Gennes et al, we use k 2 ± 
2mE/h 2 w k 2 (1 ± q/kp) 2 , with q — 2mE/h 2 k p . Using the asymptotic form 
H§ for the J v {z) ~ (\z) v /T(v + 1) we sec that the r -> limit for /(r) is 



f(r) = 



A\r v 



(D.191) 



where v = /x — 1/2, i.e. there is integer power law decay of the particle and 
hole wave functions approaching the origin. 



D.2. Asymptotic behaviour of the semiclassical wave functions 

Our two semiclassical wave functions are: 
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For the single vortex problem these become: 



with 



ui(r) 
wi(r) 

ui(r) 
vi(t) 



P(r) = 



= E/ j ( r > 



-\-ik z z+ifiO—icrz 6 



2 JS+2/i0 — ZCTz 



0. 



.4-' 



dr 



u J 0I (r) 
u o,i( r ) 



; |S^(r)+iS|(r) 



A' 



, dE(p,r) 

dr 



U 0,l( r ; ^) ] e i5 3 (r;?i) 



(D.192) 
(D.193) 



Let us investigate these in turn. 
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D.2.1. Ther^O form for f j (r) 



The spinor amplitudes in ([D.192) are given by equation (2.26). Ex- 



panding as 



0, and noting that v = eA/m = 0, we obtain 



M o,i( r ) 



0(r) 
~0(r) 



(Here we have assumed |A(r)| cx r.) We also have 
pHr) 



Pf 



±2m^/E 2 - |A(r)p 



•0(r), 



(D.194) 



from which it follows that 



1 



, 9E (p,r) 
dr 



e -iir/4 



+ 0{r% 



and since 



we have 



p^(r)dr = In 



r a ,b 



0(r 2 ), 



I Ml 



In the last two equations the phase reference point r a (r&) is chosen to be 
the classical turning point defined by p + (r a ) = (p~ (n,) = 0). We have not 
yet calculated Si (r) , but if we assume it can be neglected the appropriate 
superposition as r — > 0, call it fn{r), takes the form 



Mr) 



-iir/i / A + ( r / r w+1/2 \ 

■ l( ( : / / ::j,« / +o ^ +5,2 >- 



(M 1/2 



In particular we note the half integer power law dec ay of t he particle and 
hole wave functions in contrast to the exact result ( D.191 ). This is com- 
pletely unsatisfactory. One then hopes that S{(r) corrects this deficency. 
That this is indeed the case can be shown since for this problem S[ (r) can 
be calculated. Using — e _1 j0(r) ~pg/mr We have: 

ft a 



iSf(r) 



ipo 



2mr 2 

r(t a ,b) 



dt, 

pg dt 



r(t) 



2mr 2 dr^ 



dr d 
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and using ^5r- = ±(ipg/mr)(l + 0(r 2 )) 



iSf(r) = In 



r a ,b 



Tl/2 



Ja,b 



Tl/2 



+ 0(r 2 ), 



0(r 2 ^ 2 ) 



Then fn(r) becomes 
Mr) 



(D.195) 



which has the correct integer power law decay as required. By following 
through this analysis we have discovered the importance of the first order 
phase, S{(r). We have also checked that the semiclassical theory carried 
out correctly gives not only a quantisation rule but also the wave function 
as well. Does our effective semiclassical theory do as well? 



D.2.2. The r -> form for f j (r;ti) 



The spinor amplitudes in ( JD.193 ) are given by equation ( 3.44 ). Ex- 
panding these as r — > we obtain 



w+j(r;ft) 
u'^r-h) 



+ 0{r 3 ) 
+ 0{r 3 ) 



The r 3 rather than r correction arises since v s = —h/2mr so that p • v s = 
—hpg/2mr 2 . 
Now 



Pr(r;H) 



\ 



„ h 2 (u 2 + l/A) 
P 2 F -Pl I J - ±2m\l I E 



^) -|A(r)p, 



2mr 2 



(D.196) 



Notice the appearance of |/i=Fl/2| here rather than \fi\ in equation ( D.194 ) 
It follows from this that 



1 



r/4 



0(r 2 ), 



dr 
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and 



/ \ |PT1/2| 
\ r a,b J 

Thus the appropriate combination, call it f?(r; K), takes the form 



t/4 / ,+ (^-1/2) / / x„ \ 

M^) = ^[ A .^ ( ^L )+«*»>. 

which has the correct integer power law decay agreeing with the previous 
two solutions. In the present case the inclusion of the 1/4 in equation 
( |D.196| ) was crucial to obtain the dependence upon |/i ± 1/2 1 . Th is 1/ 4 
appears in the pf 1 due to the inclusion in the Hamiltonian, equation ( 5.90| ) , 



of h 2 terms. This confirms that our procedure, differing from that proposed 
by Littlejohn and Flynn fl7j] , is the correct one for a semiclassical theory 
of superconductors. 

We have successfully verified both versions of the semiclassical wave 
function derived in this paper. 
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